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PREFACE 


This  research  was  initiated  in  1972  during  the  author's  final  year  of 
residence  at  the  University  of  Colorado,  Boulder,  Colorado.  The  work  was 
undertaken  as  partial  fulfillment  of  the  requirement  for  a Doctor  of  Philosophy 
degree  in  Aerospace  Engineering  Sciences  at  the  University  of  Colorado. 

The  major  portion  of  the  research  was  completed  while  the  author  was  assigned 
to  the  Frank  J.  Seiler  Research  Laboratory  (AFSC) . The  computational  effort  was 
accomplished  primarily  at  the  Education  and  Research  Computer  Center  facilities 
of  the  United  States  Air  Force  Academy. 

The  perturbation  technique  used  in  this  research  was  pursued  sufficiently 
to  reveal  a more  simplified  procedure  for  rapid  recovery  of  the  coefficients 
which  must  be  calculated  in  the  derivation  of  the  perturbation  equations.  The 
original  procedure  tor  these  calculations  required  an  extensive  period  of  tedious 
effort  to  produce  the  perturbation  expressions  through  fourth  order  in  successive 
approximation.  The  revelation  of  the  simplified  procedure  made  it  possible  to 
reproduce  the  first  four  orders  of  results  in  significantly  reduced  time  as  well 
as  calculate  an  additional  five  orders  for  the  perturbation  equations. 


This  report  documents  those  results  obtained  from  the  initial  effort  to 
derive  the  first  four  orders  in  the  expansion  employing  the  original  perturbation 
computational  methods. 
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CHAPTER  I 


INTRODUCTION 


The  structure  of  a normal  shock  wave  in  a monatomic  gas  has 
been  investigated  extensively  by  theoreticians  and  experimentalists 
for  many  years.  The  shock  wave  may  be  represented  as  a one  dimen- 
sional gasdynamic  phenomenon  in  which  large  departures  from  thermo- 
dynamic equilibrium  occur  over  a few  mean  free  paths  within  the 
shock  wave.  The  one  dimensional  shock  wave  should  be  completely 
described  by  the  Boltzmann  equation,  if  the  Boltzmann  equation  is 
correct.  But,  because  of  the  nonlinearity  and  complex  nature  of  the 
collision  integral,  trustworthy  solutions  of  the  Boltzmann  equation 
have  been  obtained  only  in  cases  close  to  equilibrium,  and  certainly 
not  for  strong  shock  waves. ^ 

This  effort  is  a theoretical  investigation  of  the  higher 
order  gasdynamic  equations  for  shock  wave  structure  obtained  by 
applying  the  well  known  Chapman-Enskog  development  to  the  Boltzmann 
equation  through  fourth  order  in  spatial  derivatives.  The  motivation 


The  relative  change  in  macroscopic  variables  over  a local 
mean  free  path  becomes  appreciable  even  for  modest  upstream  Mach 
numbers  (M^) . For  example,  for  M.  = 1.5  the  relative  change  in 
macroscopic  variables  over  a local  mean  free  path  is  less  than  0.1, 
but  for  = 2.5  the  relative  change  in  macroscopic  variables  over 
a local  mean  free  path  exceeds  0.4  in  the  interior  of  the  shock 
wave.  These  estimates  are  based  on  numerical  integration  of  the 
Navier-Stokes  differential  equations. 
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of  this  work  is  to  strive  toward  a test  of  the  nonlinear  features  of 

the  Boltzmann  equation,  which  has  not  been  done  before.  The 

Chapman-Enskog  development  has  been  used  because  it  has  been  shown 

to  have  quantitative  physical  significance  beyond  the  Navier-Stokes 

2 

level  in  studies  of  high  frequency  sound  propogation. 

An  additional  motivation  for  using  the  Chapman-Enskog 
development  has  been  the  desire  to  reconcile  conflicting  opinions 
in  the  literature  about  its  utility  drawn  from  studies  of  high 
frequency  sound  waves'^  and  weak  shock  waves. ^ 

The  gasdynamic  equations  derived  and  investigated  in  this 
report  are  for  a gas  of  Maxwell  molecules  (inverse  fifth  power 
repulsion).  The  restriction  to  Maxwell  molecules  affords  consider- 
able mathematical  simplification  in  the  Boltzmann  collision 
operator.  Although  such  an  idealized  molecular  model  has  analytical 
advantage,  it  cannot  be  expected  to  lead  to  decisive  comparision 
with  experiment.  Nevertheless,  it  seems  reasonable  to  hope  that 
results  for  shock  waves  in  a gas  of  Maxwell  molecules  may  provide 
qualitative  insights  and  guides  to  future  work  with  a realistic 
potential,  for  such  correspondence  has  been  found  for  other  phenomena. 


2 

“J.  D.  Foch  Jr  and  G.  W.  Ford,  in  Studies  in  Statistical 
Mechanics , edited  by  J . de  Boer  and  G.  E.  Uhlenbeck  (North  Holland, 
Amsterdam,  1970),  Vol.  V. 

3 

C.  S.  Wang  Chang  and  G.  E.  Uhlenbeck,  in  Studies  in  Statis- 
tical Mechanics,  edited  by  J.  de  Boer  and  G.  E.  Uhlenbeck  (North 
Holland,  Amsterdam,  1S70) , Vol.  V. 

^L.  Talbot  and  F.  S.  Sherman,  NASA  Memorandum  12-14-58W 
(1959);  F.  S.  Sherman  and  L.  Talbot,  in  Proceedings  of  the  First 
International  Symposium  on  Rarefied  Gas  Dynamics,  edited  by  F.  M. 
Devienne  (New  York:  Pergamon,  1960),  161. 


The  Chapman-Enskog  expansion  applied  to  the  Boltzmann 
equation  yields  successively  the  Euler  (zeroth  order),  Navier- 
Stokes  (first  order),  Burnett  (second  order),  Super-Burnett  (third 

order),  Super-Super-Burnett  (fourth  order),  gasdynamic 

differential  equations. 

The  derivation  and  solution  of  the  Navier-Stokes  equations 
are  well  documented  in  the  literature.^  Gilbarg  and  Paolucci^ 
successfully  integrated  the  first  order  equations,  thus  providing  a 
practical  method  for  constructing  shock  profiles  from  the  differ- 
ential equations.  Standard  Runge-Kutta  numerical  integration 
methods  yield  results  for  the  variation  of  gas  properties  through 
the  shock  for  any  incident  Mach  number  (M^) . 

The  results  of  similar  efforts  applied  to  the  Burnett  differ- 
ential equations  were  reported  by  Sherman  and  Talbot, ^ who  found  that 
numerical  integration  was  successful  only  for  Mach  numbers  less  than 

g 

about  2.1.  One  of  the  early  findings  in  the  present  work  was  that 


^S . Chapman  and  T.  G.  Cowling,  The  Mathematical  Theory  of 
Non-Uniform  Gases  (Third  edition,  Cambridge  University  Press, 

1970);  W.  G.  Vincenti  and  C,  H.  Kruger  Jr.,  Introduction  to  Physical 
Gas  Dynamics  (John  Wiley  & Sons,  1965);  L.  M.  Schwartz  and  D.  F. 
Hornig,  Physics  of  Fluids  L2,  1669  (1963);  D.  Gilbarg  and  D. 
Paolucci,  Journal  of  Rational  Mechanics  and  Analysis  2^,  617  (1953). 

^Gilbarg  and  Paolucci,  Journal  of  Rational  Mechanics,  op. 

cit. 

^ Talbot  and  Sherman,  NASA  Memorandum,  op.  cit.;  Sherman 
and  Talbot,  Proceedings  of  the  First  International  Symposium,  op. 
cit. 


Q 

The  equations  used  by  Sherman  and  Talbot  contain  an  error, 
but  the  error  has  negligible  effect  on  their  results.  In  the  expres- 
sion for  p on  page  28  of  reference  2,  and  again  on  page  164  of 

reference  3,  the  coefficient  of  (u^/p) • (du/dx)^  is  40/27;  it  should 
be  8/9. 


-t  i^TjS****" 


the  Super-Burnett  and  Super-Super-Burnett  equations  are  not  amenable 

to  conventional  methods  of  numerical  integration  for  any  Mach  number. 

9 

Nevertheless,  Montgomery  has  reported  a proof  of  the  existence  of 
shock  wave  solutions  for  the  higher  order  gasdynamic  equations  for 
Mach  numbers  close  to  unity,  although  the  uniqueness  of  these  solu- 
tions is  not  yet  established. 

Recent  experimental  findings  of  Alsmeyer^  and  earlier 
findings  of  Schmidt^  provide  evidence  of  the  inadequacy  of  the 
Navier-Stokes  equations  for  describing  shock  wave  structure.  The 
discrepancy  between  the  first  order  Navier-Stokes  theory  and  experi- 
ment becomes  conspicuous  when  the  asymmetry  of  the  normalized 
density  profile  R^(z)  is  examined  as  a function  of  Mach  number  M^. 
Following  the  conventions  of  Schmidt  and  Alsmeyer,  ^p(z)  anc*  Qp  are 
defined  according  to  the  following  relationships: 


p(z)-p. 

R (z)  = — 

P P2-P1 


(1-1) 


fIR . (z)  ] d z 

0 = =&& (1-2) 

P /[ 1-R  (z)]dz 

o p 

In  equation  (1-1)  the  subscript  1 corresponds  to  conditions  far 
upstream  of  the  shock  wave,  and  the  subscript  2 corresponds  to 
conditions  far  downstream  of  the  shock  wave.  This  notation  and  the 


^J.  T.  Montgomery,  Physics  of  Fluids  18^,  148  (1975). 

^H.  Alsmeyer,  to  be  published  in  Journal  of  Fluid  Mechanics. 
^B.  Schmidt,  Journal  of  Fluid  Mechanics  3£,  361  (1969). 


choice  of  the  origin  of  the  coordinate  system  at  the  center  of  the 
density  profile  (R  =0.5)  are  illustrated  in  Figure  1. 


Figure  1 - Characteristic  variation  of  the  normalized  density 
across  a normal  shock  wave. 

The  Navier-Stokes  theoretical  result  for  based  on  a 

0 68 

realistic  potential  (viscosity  y a T ) and  the  data  from  electron 
beam  measurements  of  R^(z)  (Schmidt  1969,  Alsmeyer  1975)  are  plotted 
versus  the  upstream  Mach  number  in  Figure  2.  It  is  quite  evident 
from  Figure  2 that  the  Navier-Stokes  theory  gives  unsatisfactory 
results  even  at  low  Mach  numbers. 

The  failure  of  numerical  integration  methods  to  provide 
solutions  to  the  higher  order  gasdynamic  differential  equations 
necessitated  the  development  in  this  effort  of  an  analytic  (pertur- 
bation) method  to  obtain  the  solutions.  The  gas  properties  were 
sought  as  perturbation  expansions,  with  the  expansion  parameter  a 
suitable  function  of  the  shock  Mach  number.  These  perturbation 
methods  were  applied  to  yield  successive  approximations  for  gas 
property  variations  across  the  shock. 

The  numerical  integration  solutions  for  the  Navier-Stokes  and 
Burnett  equations  served  as  test  cases  for  the  perturbation 
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expansions.  Reasonable  agreement  between  the  results  of  the  two 
approaches  was  obtained  for  the  Navier-Stokes  and  Burnett  equations 
and  this  agreement  served  as  the  rationale  for  carrying  out  the 
perturbation  calculations  for  the  higher  order  gasdynamic  equations. 

In  addition  to  the  authors  already  mentioned,  some  additional 
•references  are  appropriate  here  to  cite  other  theoretical  and  experi- 
mental efforts  of  relevance  to  this  report.  No  attempt  will  be 
made  to  give  an  historical  account  of  the  extensive  efforts  of  the 
many  investigators  of  shock  wave  structure.  For  such  an  account, 

the  reader  is  referred  to  a stimulating  and  thorough  publication  by 
12 

Fiszdon  containing  an  extensive  bibliography  of  some  150  theoreti- 
cal and  experimental  references  on  the  subject.  In  his  paper, 

Fiszdon  classifies  the  theoretical  attempts  to  obtain  solutions  for 
shock  structure  as  a)  Continuum  Gas  Theories,  b)  Kinetic  Models  of 
the  Boltzmann  Collision  Integral,  c)  Kinetic  Models  of  the  Distri- 
bution Functions,  d)  Monte  Carlo  Methods. 

In  the  first  theoretical  class,  the  principal  method  is  the 
application  of  the  Chapman-Enskog  development  to  the  Boltzmann 

equation  to  obtain  the  constitutive  gasdynamic  relations  for  sh)ck 

13 

structure.  In  the  second  class,  the  BGK  method  is  one  of  the  best 
known  efforts  to  develop  a simplified  model  of  the  Boltzmann 

12 

W.  Fiszdon,  in  Proceedings  of  the  Ninth  International 
Symposium  on  Rarefied  Gas  Dynamics,  edited  by  E.  Becker  (Academic 
Press,  1974)  2,  B.23. 

13 

P.  L.  Bhatnagar,  E.  P.  Gross,  and  M.  Krook,  Physical 
Review  94,  511  (1954). 


14 

collision  integral.  However,  as  Alsmeyer  reports,  experimental 
results  show  clearly  the  failure  of  the  BGK  solutions  even  for  low 
Mach  numbers. 

The  bimodal  velocity  distribution  function  introduced  by 
Mott-Smith''""*  and  the  subsequent  shock  thickness  calculations  of 
Muckenfuss^  are  characteristic  of  the  kinetic  model  investigations 
in  the  ^hird  theoretical  class.  The  Mott-Smith  method  yields 
symmetric  density  profiles  at  all  Mach  numbers,  whereas  the  experi- 
mental density  profiles  are  asymmetric  for  all  Mach  numbers. 

In  recent  years,  some  of  the  most  interesting  studies  of  shock 
structure  have  centered  on  the  development  of  Monte  Carlo  techniques 
in  an  effort  to  find  a solution  of  the  complete  Boltzmann  equation. 

Hicks  and  Yen  and  their  associates  ^ have  used  Nordsieck's  Monte 
18 

Carlo  method  to  "solve"  the  Boltzmann  equation  for  the  plane  steady 


^Alsmeyer,  op.  cit. 

"^H.  M.  Mott-Sraith,  Physical  Review  <82,  885  (1951) . 

^ C.  Muckenfuss,  Physics  of  Fluids  _5,  1325  (1962). 

^ A.  Nordsieck  and  B.  L.  Hicks,  in  Proceedings  of  the  Fifth 
International  Symposium  on  Rarefied  Gas  Dynamics,  edited  by  C.  L. 
Brundin  (Academic  Press,  1967)  1^,  695;  S.  M.  Yen,  Physics  of  Fluids 
j),  1417  (1966);  B.  L.  Hicks  and  S.  M.  Yen,  Physics  of  Fluids  10,  458 
(1967);  T.  Holtz,  E.  P.  Muntz,  and  S,  M.  Yen,  Physics  of  Fluids  14 , 

545  (1971);  S.  M.  Yen,  International  Journal  of  Heat  and  Mass  Trans- 
fer  W,  1865  (1971);  B.  L.  Hicks,  S.  M.  Yen,  and  J.  Reilly,  Journal 
of  Fluid  Mechanics  53^,  85  (1972);  S.  M.  Yen,  W.  P.  Walters,  W.  Ng, 
and  J.  R.  Flood,  in  Proceedings  of  the  Eighth  International  Symposium 
on  Rarefield  Gas  Dynamics,  edited  by  K.  Karamcheti  (Academic  Press, 
1974)  1.,  137;  S.  M.  Yen,  W.  Ng,  R.  M.  Osten,  and  W.  P.  Walters, 
Coordinated  Science  Laboratory  Progress  Report,  July  1971  through 
June  1972,  August  1972;  S.  M.  Yen  and  W.  Ng,  Journal  of  Fluid  Mechanics 
65,  127  (1974). 

18 

Nordsieck  and  Hicks,  op.  cit. 
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shock  wave.  In  an  independent  effort.  Bird  has  developed  a Monte 

Carlo  direct  simulation  method  or  numerical  experiment  in  which  he 

tracks,  in  physical  space,  the  motions  and  interactions  of  the  mole- 

20 

cules  across  the  shock.  Bird  has  claimed  that  his  direct  simulation 
method  is  entirely  consistent  with  solving  the  Boltzmann  equation. 

The  manifestly  different  stochastic  techniques  called  Monte 
Carlo  methods  by  Hicks,  Yen  and  associates  and  by  Bird  yield  substan- 
tially different  quantitative  results  for  the  density  profiles  for 
both  elastic  spheres  and  Maxwell  molecules.  This  conspicuous  differ- 
ence in  the  results  for  each  of  the  molecular  models  is  clearly 

21 

illustrated  in  the  references  and  is  attributed  to  the  differences 

in  the  Monte  Carlo  method  used. 

The  Mont'  ~arlo  methods  have  been  reported  to  give  the  best 

agreement  with  experiment.  In  addition  to  the  elastic  sphere  and 

Maxwell  molecule  data  calculated  by  Hicks,  Yen  (et  al.)  and  by  Bird, 

Bird's  method  has  been  used  to  canculate  shock  profiles  for  a number 

22 

of  other  repulsive  intermolecular  potentials.  Schmidt  reported 
excellent  agreement  between  his  experimental  results  in  argon  and 

19G.  A.  Bird,  Physics  of  Fluids  6,  1518  (1963);  G.  A.  Bird, 
in  Proceedings  of  the  Fourth  International  Symposium  on  Rarefied 
Gas  Dynamics,  edited  by  J.  H.  de  Leeuw  (Academic  Press,  1965)  1., 

216;  G.  A.  Bird,  Journal  of  Fluid  Mechanics  30,  479  (1967);  G.  A. 

Bird,  in  Proceedings  of  the  Sixth  International  Symposium  on  Rarefied 
Gas  Dynamics,  edited  by  L.  Trilling  and  H.  Y.  Wachman  (Academic  Press, 
1969)  1,  301;  G.  A.  Bird,  Physics  of  Fluids  13,  1172  and  2676  (1970). 

20 

Bird,  Physics  of  Fluids  13.,  op.  cit. 

21 

W.  Fiszdon,  op.  cit.;  Yen  and  Ng,  op.  cit. 

22 

Schmidt,  op.  cit. 


Bird's  inverse  twelfth-power  force  data  at  one  Mach  number  (M^  = 8.0). 

23 

In  later  experiments,  Alsmeyer  indicates  excellent  agreement 
between  his  experimental  data  and  the  Bird  Monte  Carlo  calculations 
for  the  inverse  tenth-power  force  data  for  the  complete  range  of 
Mach  numbers  up  to  = 9.0.  Alsmeyer  also  suggests  that  it  is 
reasonable  to  expect  improved  results  from  the  Hicks  and  Yen  method 
for  more  realistic  intermolecular  potentials. 

In  a recent  effort  to  infer  the  intermolecular  potential 

24 

from  shock  structure,  Sturtevant  and  Steinhilper  determined 
intermolecular  potential  parameters  for  four  monatomic  gases  from 
accurate  shock  tube  experiments  using  Bird's  Monte  Carlo  results. 

The  results  indicate  that  the  profiles  are  generally  insensitive  to 
different  potentials  at  high  Mach  numbers.  Reliable  results  at 
low  Mach  numbers  are  not  yet  available. 

The  present  efforts  to  derive  and  solve  higher  order  gas- 
dynamic  equations  of  shock  structure  will  be  discussed  briefly  to 
summarize  the  developments  given  in  the  succeeding  chapters  of  this 
report. 

Chapter  II  includes  the  derivation  of  the  higher  order  differ- 
ential equations  for  the  density  (or  velocity)  and  temperature  as  func- 
tions of  distance  through  the  shock  wave.  The  conservation  equations 

23 

Alsmeyer,  op.  cit. 

^ B.  Sturtevant  and  E.  A.  Steinhilper,  in  Proceedings  of  the 
Eighth  International  Symposium  on  Rarefied  Gas  Dynamics,  edited  by 
K.  Karamcheti  (Academic  Press,  1974)  .1,  159. 


are  given  for  conditions  inside  the  shock  and  the  higher  order 
constitutive  relationships  are  obtained  from  the  Chapman-Enslcog 
expansion  applied  to  the  Boltzmann  equation.  The  derivation  of  the 
differential  equations  in  each  spatial  order  involves  the  deter- 
mination of  the  corresponding  heat  flux  and  viscous  stress  contri- 
butions. The  highest  order  of  the  spatial  derivatives  determined 
in  these  contributions  specify  the  order  of  the  corresponding  differ- 
ential equations. 

Chapter  III  includes  methods  to  solve  the  Navier-Stokes 
equations  of  shock  structure.  The  standard  numerical  integration 
methods  are  presented,  followed  by  the  introduction  of  the  perturba- 
tion technique  which  was  developed  to  obtain  solutions  of  the  higher 
order  differential  equations.  The  perturbation  method,  though  not 
specifically  required  to  provide  solutions  of  the  Navier-Stokes 
equations,  is  introduced  in  Chapter  III  because  the  rationale  for 
its  applicability  to  the  higher  order  developments  is  the  same  as 
that  applied  to  the  Navier-Stokes  equations. 

Chapter  IV  provides  a discussion  of  the  limitations  of 
standard  numerical  integration  methods  in  the  Burnett  and  higher 
order  developments.  The  results  of  the  higher  order  perturbation 
calculations  are  presented,  and  a pattern  is  discerned  among  the 
resulting  equations  which  reveals  a formal  similarity  among  the 
Navier-Stokes  and  higher  order  results. 

In  Chapter  V,  shock  profiles  obtained  from  the  numerical 
integration  and  perturbation  methods  are  presented  and  compared. 


r 


12 

The  utility  of  the  perturbation  method  is  evident  from  the  relative 
agreement  between  the  results  of  the  two  methods  in  the  Navier- 
Stokes  and  Burnett  (M^<2.1)  developments. 

Chapter  VI  provides  a qualitative  look  at  the  present  theo- 
retical results  for  Maxwell  molecules  in  light  of  recent  experimental 
findings.  The  potential  need  and  means  for  extending  the  pertur- 
bation method  to  higher  orders  are  discussed  as  an  outlook  for  future 


work. 


CHAPTER  II 


HIGHER  ORDERS  OF  THE  CHAPMAN-ENSKOG 
DEVELOPMENT  FOR  MAXWELL  MOLECULES 

Conservation  Equations  for  the  Steady  State 
One  Dimensional  Normal  Shock  Wave 

The  higher  order  gasdynamic  equations  for  a Maxwell  gas  will 
be  derived  from  the  Boltzmann  equation  in  this  chapter  as  a first 
step  in  the  calculation  of  the  structure  of  a normal  shock  wave. 

The  derivation  is  patterned  after  the  well  known  Chapman-Enskog 
development  for  solving  the  Boltzmann  equation  in  successive  approxi- 
mat ion . 

The  starting  point  for  the  Chapman-Enskog  development  is  the 
Boltzmann  equation  for  the  velocity  distribution  function  f(r,v,t), 

If  = ~ vi  + /v1dVdfi80(g’x)tf,fi  " ffi]>  (II-1) 

The  first  term  on  the  right  gives  the  rate  of  change  of  molecules 
of  velocity  class  v which  are  contained  in  a specific  volume  element 
of  physical  space  with  velocity  in  a specific  volume  element  of 
velocity  space  due  to  free  flight.  The  integral  on  the  right-hand 
side  provides  the  rate  of  change  of  molecules  of  velocity  class  v 
due  to  collisions  with  other  molecules.  The  integral  sum  of  all 
encounters  with  molecules  (specified  as  class  v^)  includes  the 
collisional  processes  which  either  replenish  or  deplete  the  molecular 


count  in  the  volume  elements  of  the  physical  and  velocity  space.  Ihe 
velocity  distribution  functions  f and  f^  in  the  collision  integral 
are  representative  of  conditions  both  prior  to  (unprimed)  and  after 
(primed)  collision. 

In  the  collision  integral  g is  |v  - v^|,  and  a( g,x)  is  the 

differential  cross  section  for  a.i  encounter  in  which  the  relative 

velocity  vector  is  turned  or  scattered  through  an  angle  x into  the 

element  of  solid  angle  dft.  For  a gas  consisting  of  atoms  of  mass 

m interacting  with  an  intermolecular  potential  4>(r)  = the 

r 

expression  for  the  cross  section  is  written  as 


go(g,x)  = (~)2  F(X)* 


(II-2) 


Such  a (fictitious)  gas  was  considered  by  Maxwell,  and  is  referred 
to  as  the  Maxwell  model. 

To  develop  the  Chapman-Enskog  expansion  for  departures  from 
the  equilibrium  condition,  the  velocity  distribution  function  is 


sought  in  the  form 


f = f(0)  + f(1)  + f(2)  + 


( 11—3) 


where  f^  is  the  local  Maxwell  distribution,  f^  is  proportional 

(2) 

to  spatial  derivatives,  f is  proportional  to  squares  of  spatial 
derivatives  and  second  derivatives,  etc.  The  local  equilibrium 
distribution  f^  is  given  by 


:<0) 


n( 


2trkT 


) exp{- 


m 

2kT 


<J  - S)2) 


( 1 1— 4 ) 


where  m is  the  molecular  mass,  n is  the  local  molecular  number  den- 
sity, T is  the  local  gas  temperature,  and  k is  Boltzmann's  constant. 


The  expression  (v  - u)  is  called  the  peculiar  or  thermal 

. ->■ 
velocity  of  the  gas,  is  denoted  by  V,  and  is  equal  to  the  difference 


between  the  molecular  velocity  (v)  an^  the  mean  flow  velocity  of  the 

m2 

molecules  (u) . The  combination  has  dimensions  of  (1/velocity) 

and  will  be  employed  frequently  m nondimensionalizing  velocity  terms 
throughout  this  work. 

Using  the  Boltzmann  equation,  as  presented  above,  and  taking 

1 2 

moments  with  respect  to  v of  1,  v,  and  y v yields  the  one-dimensional 
conservation  equations  for  the  normal  shock  wave.  The  mean  gas  flow 
direction  is  taken  along  the  Z-coordinate.  Detailed  treatments  of 
the  moment  calculations  are  plentiful  in  the  literature'*'  and  will  not 
be  repeated  in  this  report.  The  results  are  given  below. 

MASS: 


3n  , 3n  3u  A 

31  + u 37  + n 37  = 0 


(II-5) 


MOMENTUM: 

3u 


__  3u  . _L  3P 
3t  u3z  nm  3z 


(II— 6) 


ENERGY : 


2 nk  (9I  + u ii)  + 2a  + P = o 

2 nK  '3t  3z;  3z  3z 


(II— 7 ) 


J.  D.  Foch  Jr.  and  G.  W.  Ford,  in  Studies  in  Statistical 
Mechanics , edited  by  J.  de  Boer  and  G.  E.  Uhlenbeck  (North  Holland, 
Amsterdam,  1970),  Vol.  V;  S.  Champman  and  T.  G.  Cowling,  The 
Mathematical  Theory  of  Non-Uniform  Gases  (Third  edition,  Cambridge 
University  Press,  1970);  W.  G.  Vincent!  and  C.  H.  Kruger  Jr., 
Introduction  to  Physical  Gas  Dynamics  (John  Wiley  & Sons,  1965). 
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In  these  equations 

P = P (z,t)  = JdVfm  V2  (11-8) 

and 

q = q (z,t)  S JdVf  | V2Vz  (II-9) 

where  P is  the  pressure  minus  the  viscous  stress  and  q is  the  heat 
flux. 


The  steady  state  equations  for  a normal  shock  wave  are  easily 
found  from  the  above  by  neglecting  the  time  derivatives  and  noting 
that  the  momentum  equation  can  be  used  in  the  energy  equation  to  form 
an  exact  energy  differential  equation.  Since,  in  steady  state,  all 
three  differential  equations  are  exact,  each  equation  can  be  inte- 
grated once  to  give  the  results  below. 

MASS: 


pu  = 

MOMENTUM: 


2 2 
pu  + P = p^u^  + n^kT^ 


ENERGY : 


[ 3 k u2) 

'5  kT  , Ui) 

[im  T + Tjpu  + uP  + q - 

--  i T 4.  

2 m i 2 

piUi 


(II-5a) 


(II-6a) 


(II-7a) 


In  these  equations  p is  the  mass  density.  The  quantities  on  the 
right-hand  sides  of  the  equations  are  integration  constants,  evalu- 
ated either  upstream  (i  = 1)  or  downstream  (i  = 2),  which  represent 
the  constant  fluxes  of  mass,  momentum,  and  energy.  Equations 
(II-5a),  (II- 6a),  and  (II- 7a)  describe  the  shock  wave  illustrated  in 


¥ 


it 


Figure  1. 


Alternate  Forms  of  the  Viscous  Stress  and  Heat 
Flux  for  the  Higher  Order  Developments 
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The  uniformity  expansion  for  the  velocity  distribution 

function  induces  uniformity  expansions  for  the  heat  flux  q^  and 

stress  P..,  which  reduce  in  the  one-dimensional  case  to  q and  P , 
ij’  Hz  zz 

or  simply  q and  P.  Following  methods  employed  in  the  Chapman- 
Enskog  solution  of  the  Bolt2mann  equation,  the  velocity  distribution 
function  is  further  defined  in  the  following  expanded  series  form: 


; (0) 


i + (j/1)  + <(>(2)  + <J>(3)  + 

r (0) 


(11-10) 


In  zeroth  order,  f = f , which  implies  q = 0 and  P = nkT,  and  the 
conservation  equations  reduce  to  the  equations  of  equilibrium  flow. 
The  inclusion  of  the  higher  order  contributions  to  f leads  to  the 
Navier-Stokes , Burnett,  and  higher  order  macroscopic  equations  of 
gas  dynamics. 

Again  following  Chapman  and  Enskog,  the  following  require- 
ments will  be  imposed  in  the  development: 


Jdvf(n)  =0,  n>l 

(11-11) 

jdvl^n\z  = 0,  n>_l 

(H-12) 

Jdvf^v^  = 0,  n>l. 

(11-13) 

These  requirements  are  consistent  with  the  equilibrium  solution  but 
are  not  necessitated  by  it.  The  requirements  will  introduce  some 


tNote:  For  clarity  in  expressing  the  results  of  the  moment 
calculations  throughout  the  remainder  of  this  chapter,  the  arguments 
(v)  and  (v^)  for  the  velocity  distribution  functions  have  been 
omitted  in  the  equations. 


very  effective  simplifications  in  the  host  of  moment  calculations 
which  are  required  to  determine  the  gasdynamic  forms  of  P and  q. 

It  is  convenient  to  express  the  integral  requirements  (11-11) , 
(11-12),  and  (11-13)  in  terms  of  the  peculiar  velocity  V: 


Jdvf(n)  = 0,  n>l 

(II-lla) 

/dvf(r)V  = 0,  n>l 

(II-12a) 

Jdfav'n)V  = 0,  n>l. 

(II--13a) 

Introducing  the  expanded  form  for  the  velocity  distribution 
function  in  ( 1 1—8 ) and  (II-9) , the  latter  expressions  give 


P = nkT  + P^  + P^  + P^  + 


(1)  . (2)  , (3) 

= q + q + q 


q = q 


(H-14) 

(11-15) 


where 


P(n)  = JdVf(n)mV2,  n>l , (11-16) 

q(n)  = /dVf(n)  f V2Vz,  n>l.  (II-17) 


With  the  aid  of  (II-12a)  and  (II-13a)  the  expressions  for  (11-16)  and 
(11-17)  may  also  be  written  in  the  form 


P'“'  = - JdVf 


,(n) 

. (n) 


(n) 


r3  2 

— V - — 
2 z 2 


- /dv£(n>  f li 


■lire  El  “1  y _ \]^\J 

2 2 Z 2 


(11-18) 


(11-19) 


where  the  equality  of  dv  and  dV  in  velocity  space  has  been  employed. 

The  velocity  forms  in  (11-18)  and  (11-19) , nondimensionalized  by 

,1 

yr 

, are  eigenfunctions  associated  with  the  linearized 


factors  of 


m 


2kT 


collision  operator  for  Maxwell  molecules.  These  and  other 
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eigenfuctions  (and  their  eigenvalues)  are  tabulated  in  Appendix  A. 
\ + 

V and  the  notation  of  Appendix  A leads  to 


*>■ 

Using  c = 
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2kT 
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Spatially  Ordered  Equations  Derived  from  the 
Expanded  Velocity  Distribution  Function 
Applied  to  the  Boltzmann  Equation 

The  developments  which  follow  in  the  remaining  sections  of 

this  chapter  involve  solution  of  the  Boltzmann  equation  using  the 

expanded  form  of  the  velocity  distribution  function  through 

f^  = but  only  to  the  extent  necessary  to  determine  P 

(4)  (4) 

and  q through  P and  q 

In  applying  the  full  Chapman-Enskog  method  to  the  solution 
of  the  Boltzmann  equation,  the  usual  procedure  is  to  substitute 
expanded  forms  of  the  velocity  distribution  (i.e.,  f = f^, 
f = f^  + f^,  • * * etc.)  into  both  sides  of  the  equation  and, 
in  successive  approximation,  solve  for  the  Subsequently,  the 

<j>  are  used  in  the  integrals  which  define  P^  and  q^n\ 

The  present  effort  has  been  designed  to  circumvent  complete 
determination  of  the  since  the  required  expressions  for  the 

viscous  stress  and  heat  flux  contributions  may  be  obtained  by  less 
cumbersome  means. 

In  the  present  effort,  as  in  the  Chapman-Enskog  method,  an 
essential  point  is  the  elimination  of  substantial  derivatives  of 


density,  flow  velocity,  and  temperature  in  favor  of  spatial  deriva- 
tives using  conservation  equations  (II— 5) , (II-6) , and  (II— 7 ) . It 
is  this  point,  as  will  become  apparent  subsequently,  that  leads  to 
dependence  of  the  f^  on  spatial  derivatives. 

By  suitably  expanding  the  expressions  involving  f and  f^  in 
the  collision  integral,  contributions  involving  the  may  be 

separated  according  to  the  order  of  their  spatial  derivatives. 

Once  the  ordering  exercise  is  complete  for  each  term  in  the 
Boltzmann  equation,  the  ordered  contributions  on  the  left  and  right- 
hand  sides  are  equated,  giving  a series  of  equations  based  on  the 
order  of  the  spatial  derivatives.  After  some  preliminary  rearrange- 
ments, the  resulting  equations  are  in  suitable  format  for  calculating 
moments  involving  specific  <j>(v),  such  as  those  occurring  in  the 
definitions  of  P and  q.  The  "selected"  <j)(v)  consist  of  only  those 

velocity  functions  which  are  indispensable  for  calculating  the 
a (n) 

and  q 

As  will  be  seen,  the  working  equations  employed  in  this 
method  have  considerable  utility  in  that  they  allow  the  calculation 
of  P^  and  q^  from  lower  order  moments  involving  only  a few 
eigenfunctions  such  as  and  H^Cc).  In  most  cases,  the  lower 

order  moments  needed  in  a given  order  are  already  known  from  the 
previous  orders.  However,  each  new  order  does  introduce  a few  new 
moments  which  must  be  evaluated.  Even  so,  the  explicit  forms  of  the 
^ need  not  be  known  for  any  of  the  higher  order  work  in  deter- 
mining the  heat  flux  and  viscous  stress  contributions.  In  fact. 


2jl 

once  the  lower  order  moments  have  been  calculated,  the  remaining 

effort  requires  only  a very  careful  algebraic  exercise. 

To  carry  out  the  specific  calculations  as  outlined  above, 

the  one-dimensional  form  of  the  Boltzmann  equation  is  written  as 

1 

Tt  + \ If  ■ F(x)  S'£i  - t£j]  - (II-22) 

where,  to  repeat  the  differential  scattering  cross  section  for 
Maxwell  molecules 

1 

o(g,X)  - ~ ^ • (II-23) 

The  time  derivative  of  the  velocity  distribution  becomes 


«.t(0>£j7+t(l>+*<2>+*<3>  + ..„] 

+ (T  + *(1>  + *(2)  + ♦«>  + .. .011^ 

and  the  streaming  term  can  be  expressed  as 


(11-24) 


3f  t x „ m=(0)  3 
v -5—  = (u  + V )f 
z dz  Z dz 


1 + + ^ + *<3>  + 


+ (u  + V ) ML  + <J>(1)  + (J)(2)  + + 


9z 

2 (11-25) 


With  <t>  of  order  n in  spatial  derivatives,  the  left-hand 
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where  — represents  the  kth  order  part  (in  the  sense  of  spatial 
derivatives)  of  the  substantial  derivative  of  the  quantity  differen- 
tiated. Similarly,  { means  all  terms  within  the  enclosure  that 
are  of  order  k. 

Substitution  of  f and  f^  in  the  collision  term  results  in  a 
similar  ordering  according  to  the  order  of  the  spatial  derivatives 
of  the  Since,  in  local  equilibrium,  the  products  of  the 

distribution  functions  before  and  after  collision  are  equal  as  a 
consequence  of  conservation  of  energy  and  momentum  (i.e., 
j(0)  j-(O)  _ f(0)f(0)}>  the  spatially  ordered  form  of  the  collision 

integral  may  be  arranged  as 


/ dvjdfi^j  P(x)[f'f1-fEl]. 


/ dv  Jdfi 
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2k 
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l l 
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Equating  the  spatially  ordered  expressions  in  (11-26)  and  (11-27) 
leads  to  the  basic  equations  for  determining  and  q^n\ 


Determination  of  P ^ and  q^ 


The  first-order  calculations  required  to  obtain  P^  and 
require  the  evaluation  of  the  derivatives  and  the  associated 
integrals  which  arise  from  the  first-order  terms  in  (11-26)  and 


(11-27). 

In  preparation  for  the  evaluation  of  the  substantial  and 
spatial  derivative  expressions  which  will  be  used  throughout  the 
remainder  of  this  chapter,  the  following  preliminary  relationships 
will  be  helpful. 


Starting  with  the  equation  for  the  equilibrium  distribution 


,(0)  m J/4  ) m 2 2 , 
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the  substantial  derivative  — may  be  conveniently  obtained  from 
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with 
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The  resulting  substantial  derivative  may  be  expressed  in  the  form 
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Converting  to  dimensionless  velocity  forms  and  recognizing  from  the 
tables  for  the  4*  ^(c)  in  Appendix  A that 

VAn(c)  - 1,  (Il-32a) 


00' 

= ft  ■ c?' 


^01 (c)  = V 


(II-32b) 

(II-32c) 


the  substantial  derivative  may  be  further  expressed  as 
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Similar  operations  on  the  remaining  spatial  derivative  term 
of  equation  (11-26)  yields 
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The  above  result  requires  some  "construction"  of  the  terms  ^r^(c) 
from  the  dimensionless  velocity  forms. 

To  assist  in  examining  the  order  of  the  contributions  in 
the  conservation  equations  as  applied  to  the  substantial  deriva- 
tives in  (11-33),  the  conservation  equations  are  written  below  so 
as  to  isolate  the  substantial  derivatives  on  one  side  of  the 
equation.  For  later  reference,  the  spatial  order  of  each  term  is 
entered  above  the  term,  taking  into  account  the  final  order  of  the 
term  after  carrying  out  the  differentiations  and  multiplications 
indicated . 
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(1)  , (2)  , (3)  , (4)  , 
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The  first-order  part  of  (11-26)  is  constructed  from  the 
first-order  part  of  (11-33)  and  this  involves  contributions  from  all 
three  of  the  above  equations.  Summing  the  contributions  from  (11-33) 
and  (11-34)  yields 
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The  remaining  work  to  determine  the  explicit  forms  of  Pv 
and  involves  the  evaluation  of  the  first-order  part  of  the 
collision  integral  given  in  (11-27),  which  is  written  below  as 


Jv  dvjdfiff)172?  (X)  f (0)  (v)f^0)  (vx)  A <f>(1)  , 

^ > • 


(11-39) 


where 


- *a>-  + _ +a>  _ *&>, 


dVl  = dVl* 


■ (0) 


and  f (v)  is  independent  of  the  integration  over  v^. 

Expression  (11-39)  is  f^(v)  times  the  linearized  collision 
operator  J which,  when  operating  on  the  various  4*r^(c)  introduced 
earlier,  yields 


d *rl<c>  = £r£  ¥r£(c> 


Here,  J has  been  defined  by 


t a/2  , r -i 

<jij=  /Vidy<i!)(fj  F(X)f‘0)4|fr!t(j)  . 

Appendix  A provides  the  following  results: 
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(11-42) 


In  the  above  k^  is  a pure  number  determined  from 

A2  =/o  dx  sinXd-cos^)  F(x) 

and  x is  the  deflection  angle  associated  with  binary  collisions. 
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A comparison  of  (11-41)  and  (11-42)  with  (11-38)  shows  that 
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where  e ^ = -ttA2  and  = - y ttA^  are  the  eigenvalues  corresponding 
to  ^11^  and  ^02  ^ ’ 

Knowing  the  first-order  contributions  to  the  viscous 

stress  and  heat  flux  may  also  be  determined  from  the  definitions 
(11-20)  and  (11-21).  From  the  orthogonality  of  the  V ^(c)  it  follows 


that 
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where  the  integral  is  of  a type  tabulated  in  Appendix  A. 
Introducing  the  abbreviation 
,1/2 


U - f kT 


J K 


the  first-order  viscous  stress  term  may  be  written  as 
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with  p recognizable  as  the  gas  viscosity. 

A similar  calculation  gives  the  first-order  heat  flux  term 
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(11-47) 
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with  the  thermal  conductivity  K defined  by 


4 m 


(11-49) 


(2)  (2) 
Determination  of  Pv  and  qv  ' 


The  integral  equation  for  <{>  is  composed  of  the  second- 

order  parts  of  (11-26)  and  (11-27).  Separating  the  collision 

integrals  which  involve  <f/2^  from  those  which  involve  <j>^  the 

(2)  (2) 

expressions  for  P and  q can  be  determined  from 
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In  order  to  determine  P'  without  solving  explicitly  for 
<P  2^,  multiply  (11-50)  by  mV^ , integrate  over  all  v,  and  use  the 
H-Theorem  transformations  on  the  collision  integrals.  These 


operations  yield 
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where,  for  any  'Kv), 

A [¥]  = [4"  + ^ - V - . 


(II-51) 


Using  the  explicit  results  for  the  relationship  between  V', 
and  V,  (from  the  dynamics  of  binary  collisions),  integration 
over  the  solid  angle  in  the  collision  integrals  of  (11-51)  may  be 
effected  with  the  result  that 
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Several  of  the  terms  of  the  right-hand  side  of  (11-52)  do  not 
survive  the  subsequent  integration  over  v and  either  because 
they  occur  as  odd  functions  of  velocity  times  the  even  function 
f^,  or  because  of  the  requirements  (II— 11) , (XI— 12)  , or  11-13). 
Thus,  the  collision  integral  on  the  left-hand  side  of  (11-51) 


reduces  to 
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The  collision  integral  on  the  right-hand  side  of  (11-51)  vanishes. 

The  remaining  integrals  in  (11-51)  containing  the  substantial 
and  spatial  derivatives  can  be  rewritten  by  noting  that  the  inte- 
gration over  velocity  space  and  the  differentiation  with  respect  to 
time  and  physical  space  are  independent,  and  hence  commute.  Thus, 
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Recalling  that  V = v - u(z,t),  which  implies  — = - r— , 
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the  preceding  results  may  be  combined  to  give 
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Since  f^  and  <J>^  are  known  explicitly,  the  remaining 
integrals  in  (11-56)  may  now  be  evaluated  explicitly.  Doing  so 
(and  then  performing  the  differentiations)  gives 
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with  p the  pressure,  p = nkT. 

(2) 

The  q calculations  proceed  in  exactly  the  same  manner 

1 2 

except  (11-50)  is  multiplied  throughout  by  -^mV  and  integrated 
over  all  v.  The  intermediate  result  from  this  operation  is 
ttS 
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Performing  the  integrations  and  extracting  the  second-order 

(2) 

contributions,  the  final  result  for  q is  found  to  be 
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(2)  (2) 

The  derivatives  of  the  number  density  n in  P and  q may 

be  eliminated  in  favor  of  derivatives  of  the  flow  velocity  u by 

using  the  integrated  continuity  equation  (nu  - constant). 

(2)  (2) 

The  results  for  P and  q , when  substituted  into  the 
general  conservation  equations  ( 1 1—6)  and  (1 1—7 ) yield  a pair  of 


(11-59) 
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coupled,  nonlinear  second-order  differential  equations  for  the 
velocity  and  temperature  within  a normal  shock  wave  in  a Maxwell 


(3)  (3) 

Determination  of  P and  q 


(3)  (3) 

The  third  order  calculation  to  determine  P and  q 

follows  the  same  pattern  as  employed  in  the  previous  order.  For 

third-order,  (11-26)  and  (11-27)  give 
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Multiplying  (11-60)  throuf.iout  by  mVz  and  integrating  over 
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all  v,  the  intermediate  result  tor  P is  obtained  as 
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The  moment  of  ’^(c)  with  respect  to  <(>^  may  be  determined 

(2) 

from  the  integral  equation  (11-50)  for  4>  , without  solving 

(ON 

explicitly  for  {}>  , in  much  the  same  way  as  the  moment  of 
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with  respect  to  <J>  ( « P ) was  determined. 

Starting  from  the  identity  (see  Appendix  A) 
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Now  multiply  (11-50)  by  H^fc),  integrate  over  all  v,  and  use  the 
H-Theorem  transformations  on  the  collision  integrals.  Proceeding  as 
in  the  preceding  section  leads  to  an  alternative  form  for  the  right- 
hand  side  of  (11-63)  or,  expressed  otherwise,  a trap  for  the  moment 
of  with  respect  to  <j>  : 
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Using  (11-64)  in  (11-61)  and  performing  the  necessary  differ- 
(3) 

entiations,  P is  found  to  be 
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In  obtaining  the  above  result  for  P , substantial  deriva- 
tives  of  the  lower  order  P^  and  q^  in  (11-61)  were  determined 


using  the  identities 
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The  third-order  contribution  to  the  heat  flux,  q , is 

1 2 

obtained  by  multiplying  (11-60)  throughout  by  v?  and  integrating 

-*■  (3) 

that  result  over  all  v.  The  intermediate  form  for  q turns  out  to 

be 
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Although  the  Pv  calculations  have  already  provided  the  result  for 

*f 

the  second-order  moment  of  y^(c)  > the  additional  second-order  moments 
of  the  eigenfunctions  Y (e)  anc*  ^12 must  sti^  be  evaluated. 
However,  methods  used  to  evaluate  the  ^03  (£)  result  may  be  used  to 
evaluate  the  additional  moments  occurring  in  (11-67).  The  calcula- 
tions  of  the  second  moments  of  (c)  and  ^32^  £ive 
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Performing  the  derivatives  and  eliminating  derivatives  of  n 
in  favor  of  derivatives  of  u yields 
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The  substitution  of  the  viscous  stress  anu  heat  flux 
contributions  through  third-order  into  the  general  conservation  equa- 
tions (II-6)  and  (11-7)  yields  a pair  of  coupled,  nonlinear  third- 
order  differential  equations  which  govern  the  spatial  dependence  of 
the  velocity  and  temperature  in  a normal  shock  wave.  The  third- 
order  results  are  known  as  the  Super-Burnett  differential  equations 
for  the  structure  of  a normal  shock  wave  in  a Maxwell  gas. 


(4)  (4) 

Determination  of  P and  q 


The  fourth-order  contributions  of  the  viscous  stress  and  heat 
flux  are  found  by  applying  the  same  moment  methods  to  the  integral 
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equation  for  <£'  ' obtained  from  (11-26)  and  11-27) 
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The  moments  with  respect  to  mV  and  -xmV  V yield  the 
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(4)  (4) 

following  intermediate  results  for  P and  q respectively: 


-U ,N1/2n  P<4>  .i|2( 

2 2{m ) 15  3z 


>A  fp(3)  (2)  (1)  7 (3)  3u 

+ Dt  [P  + P + P j + 3 P 3z 

r , ,1/9  2 

3 nm  2kT  1 2f c . (3)w 

+ 3l  "3/2  — ViCS  * T03<c) 

. TT  J 

m L J 


vP(3>  +^p(2)  vp(i)  + |nkT 

Dt  Dt  . Dt  Dt  2 


. 16  3u  (3)  7 _3_  kT  p(3) 

5 3z  ^ 2 3z  m 


(11-71) 


/dciC  d»(3)  4-  £ 


+ f^nkT 


L 3 J 4 kV  1 r+-c2.0)  fur  fts  w 

3z  | 3 n m ^3/2  Jdce  * [^20(c)  “ ^12(c 


(11-72) 


The  evaluation  of  the  third-order  moments  of  ^(c),  4'20(c)’ 

and  4^2  (c)  proceeds  as  in  the  previous  sections.  The  calculations 

(4)  (4) 

of  the  nonlinear  expressions  for  P and  q were  carried  out  to 
the  point  where  the  results  were  known  in  terms  of  substantial 
derivatives,  spatial  derivatives  and  the  lower  order  viscous  stress 
and  heat  flux  contributions.  For  reasons  which  will  become  apparent 


in  the  next  chapter,  it  was  not  necessary  to  perform  all  the 


differentiations  and  substitutions  to  reduce  these  very  involved 

expressions  to  the  spatial  derivative  formats  as  was  done  in  the 

(4)  (4) 

previous  developments.  In  fact,  only  the  terms  in  P and  q 
which  are  proportional  to  the  fourth-order  spatial  derivatives 
3^n  34u  3^T 

— 27  > — T > or  — T will  be  needed  later  in  this  report.  Neglecting 
3z4  3z^  3z^ 

all  other  terms  leads  to  the  linearized  results 
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If  the  complete  nonlinear  results  for  P and  q were 
included  with  the  lower  order  viscous  stress  and  heat  flux  contri- 
butions, the  general  conservation  equations  (IT-6)  and  (II-7)  would 
yield  two  coupled,  nonlinear  fourth-order  differential  equations 
for  the  temperature  and  velocity  in  a normal  shock  wave.  These 


fourth-order  results  are  known  as  the  Super-Super-Burnett  equations. 


CHAPTER  III 


METHODS  TO  SOLVE  THE  NAVIER-STOKES  EQUATIONS 
Numerical  Integration 

It  is  convenient  to  rearrange  the  general  conservation 
equations  (II-5a) , (II-6a) , and  (II-7a)  in  two  ways.  First,  use  the 
continuity  equation  (II-5a)  to  eliminate  p in  favor  of  u in  the 
momentum  and  energy  equations;  then,  use  the  momentum  equation  to 
eliminate  P from  the  energy  equation.  These  rearrangements  give 


3^-^-  lq]  + | (T-l)  - | mJ(0-1)2  + (o-l)  = 0,  (III-1) 

and 

cw  2 

5M^  0 52 

[P]  + j Ma)(u-1)  + (T-o)  = 0,  (III— 2) 
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where  the  dimensionless  velocity  (u)  and  the  dimensionless  temper- 
ature (t)  are  defined  by 
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with  i = 1 corresponding  to  conditions  far  upstream  and  i = 2 
corresponding  to  conditions  far  downstream.  The  quantity  is 


the  Mach  number, 
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i a.  ’ 
i 


( III— 5) 


in  which  is  the  speed  of  sound  determined  from 
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In  the  Navier-Stokes  order  of  fluid  mechanics,  the 
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expressions  for  P and  q are 
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For  Maxwelx  molecules,  u and  K are  given  by 
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In  due  course,  following  the  convention  of  other  authors, 
distance  in  the  direction  of  the  mean  flow  will  be  measured  in 
multiples  of  the  upstream  Maxwellian  mean  free  path  A^,  defined  by 
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Anticipating  the  eventual  introduction  of  A^,  the  Navier-Stokes 
equations  for  a Maxwell  gas  may  be  written  in  the  form 
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Gilbarg  and  Paolucci  provided  an  extensive  analysis  of  the 
solution  of  (II I— 12 ) and  (III-13) . Everything  depends  on  the  mathe- 
matical nature  of  the  singular  points  of  (III-12)  and  (III-13)  in  the 
auxilliary  O-T  phase  plane.  The  physically  significant  solution  of 
(111-12)  and  (II 1—13 ) may  be  represented  by  a trajectory  which  con- 
nects these  two  points.  Such  a solution  exists  and  may  be  obtained 
by  numerical  integration,  but  the  numerical  integration  must  be 
initiated  near  the  downstream  singularity  and  proceed  upstream. 
Numerical  integration  from  upstream  toward  downstream  will  necet  - 
sarily  fail. 

This  requirement  for  successful  numerical  integration  may 
be  understood  by  examining  all  the  mathematical  solutions  including 
the  nonphysical  ones  very  near  the  singular  points.  To  accomplish 
this,  linearize  the  Navier-f tokes  equations,  and  put 


T = T.  He 

l 


K.Z 

1 


u . + u e 

l 


Ki* 


(II 1—14 ) 
(II 1—15 ) 


or 


T = 1 + T e 


Ki2 


(111-16) 


42 


Ps  . C. 

O = 1 + 0 e . (III-17) 

The  characteristic  or  secular  equation  for  may  be  written  in 
the  form 


20  L2  + (9-23M2)  L,  + 6 M2  (M2-l)  = C,  (III-18) 
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The  roots  of  this  characteristic  polynomial  are  given 
graphically  in  Appendix  B along  with  a table  of  the  physically  signi- 
ficant roots  used  in  the  subsequent  calculations  of  the  gas  property 
profiles.  An  examination  of  the  characteristic  roots  shows  for  all 
Mach  numbers  that  upstream,  both  roots  are  real  and  both  are  positive, 
whereas  downstream  the  roots  are  rea1  but  of  opposite  sign.  There- 
fore, near  the  upstream  singularity,  the  most  general  solution  of 
the  linearized  Navier-Stokes  equations  involves  a linear  combination 
of  exponentials  with  the  two  positive  values  of  These  exponen- 

tials tend  to  zero  as  z-+  - giving  u-ni^  and  T-*T^.  However,  near 
the  downstream  singularity,  the  two  values  of  have  opposite  sign 
and  the  most  general  solution  diverges  as  z->  + There  is  only 
one  trajectory  for  which  u-m^  and  T-H^  as  z-»-  + 00 . 

Thus,  it  is  the  signs  of  the  roots  of  the  characteristic 
equation  (III-18)  that  determine  whether  the  Navier-Stokes  equations 
are  amenable  to  numerical  integration,  and,  if  so,  in  which  direc- 
tion through  the  shock  wave.  The  means  used  to  attain  this  insight 
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are  extremely  important,  for  they  illuminate  difficulties  encoun- 
tered later  in  attempting  numerically  to  integrate  the  higher  order 
gasdynamic  equations. 

The  results  of  numerical  integration  of  the  Navier-Stokes 
equations  will  be  discussed  later  in  this  report  in  Chapter  V. 

Perturbation  Expansion 

It  will  be  seen  shortly,  that  some  difficulties  arise  whan 
numerical  integration  is  attempted  in  the  higher  order  gasdynamic 
equations.  Thus,  at  this  point  in  the  thesis,  it  is  convenient  to 
introduce  the  development  of  a perturbation  expansion  method  which 
can  be  used  to  integrate  the  higher  order  equations  where  the 
numerical  integration  fails.  It  is  appropriate  to  discuss  the 
development  of  the  perturbation  expansion  now,  in  conjunction  with 
the  Navier-Stokes  equations  for  it  will  be  assumed  that  such  a 
method  will  be  trustworthy  if  it  gives  good  agraement  with  the 
results  of  numerical  integration  of  the  Navier-Stokes  equations. 

A successful  perturbation  expansion  presumably  must  be  based 
on  a dimensionless  parameter  that  is  always  less  than  of  order  unity 
(even  for  M^>>1).  In  addition,  the  perturbation  expansion  should 
not  lead  to  approximate  results  that  violate  the  Rankine-Hugoniot 
relations  (which  connect  upstream  and  downstream  values).  Finally, 
the  perturbation  expansion  must  somehow  do  justice  to  the  enormous 
variation  in  shock  wave  thickness  between  the  extremes  of  M.^1 


and  M^>>1. 
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In  consideration  of  the  above,  the  choice  of  the  expansion  parameter 
2 

is  e = M2  - 1,  and,  from  this  definition,  it  is  noted  that  £ varies 
monotonically  from  0 to  -0.8  as  increases  from  1 to  00 . 

In  order  to  avoid  violations  of  the  Rankine-Hugoniot  relations 
induced  by  the  perturbation  expansion,  a normalized  density  R^  and 
normalized  temperature  R^  are  regarded  as  basic  dependent  variables, 
where 


P(z)-P, 


Ro(z)  = 
p p2  M1 

T(z)-T 
R (z)  = — — — 

2"  1 


(in-21) 


(III-22) 


These  normalized  quantities  may  be  expected  to  vary  from  0 to  1, 
independent  of  M^.  Since  the  boundary  conditions  which  these  vari- 
ables must  satisfy  are  independent  of  M^,  it  is  unlikely  that  approx- 
imations to  Rp  and  Rt  will  fail  to  have  the  correct  asymptotic 
values.  From  (III-21)  and  (III-22)  and  the  Rankine-Hugoniot  rela- 
tions it  follows  that 
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Note  for  future  reference  that  one  could  have  introduced  a 


normalized  density  and  normalized  temperature  R^,  defined  by 
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The  pronounced  dependence  of  shock  wave  thickness  on  Mach 
number  suggests  the  introduction  of  a Mach  number  dependent  length 
scale  to  "tame"  the  differential  equations.  The  idea  is  to  regard 
Rp  and  R^  as  functions  of 

Z = < (M) *2  (III-29) 

with  k(M)  chosen  so  that  the  thickness  of  the  shock  wave  (in  units 

of  O does  not  show  a pronounced  dependence  on  Mach  number. 

How  does  this  lead  to  a choice  for  k(M)?  It  is  known  that 

<^z 

far  upstream  R^  and  R^  are  proportional  to  e , with  the  physically 
significant  choice  for  k (positive,  tending  to  zero  as  M^-PL), 


plal  1 


K1  UjM^  40 


Tr!  \ (23  Mj-9)  - 


(23  mJ-9)2  - 480  M2(M2-1) 


1/2 
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Far  downstream,  however,  R^  and  R^  are  proportional  to  e 1 , with 
the  physically  significant  choice  for  (negative,  tending  to  zero 
as  M^->1) , 


'wy’toNf' 
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K = P2a2  1 
2 y2M2  40 


~|(23M2-9)  - (23  M2-9)2-480  M2(M2-1)  1/2  | . (III-31) 


For  M~l,  the  two  reciprocal  length  scales  (III-30)  and  (III-31)  are 

2 1 
both  proportional  to  (M^-l),  but  for  M^»l,  and  i<2~—  * 

Now,  it  will  be  shown  that  in  the  interior  of  the  shock  wave, 

for  M >>1,  K~r7~.  To  do  this,  start  from  the  full  nonlinear  Navier- 
1 M1 
Stokes  equations, 
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set  i = 1 and  think  of  u and  T expressed  in  terms  of  R^  and  R^. 
Next,  assume  the  existence  of  an  "interior"  of  the  shock  wave  where 
Rp  and  Rt  are  both  of  order  0.5.  It  then  follows  from  ( III— 32)  and 
(III— 33)  that  within  such  a region,  for  M^>>1, 


dz  M, 


T JL 
dz  M, 


In  other  words,  K~;—  in  the  interior  of  the  shock  wave  for  M,>>1. 

M1  1 
Hence,  relying  on  the  foregoing  analysis,  the  Mach  number 

dependent  length  scale  for  the  perturbation  expansion  will  be  based 
on  (III-31).  This  choice  may  be  expected  to  be  appropriate  through- 
out the  shock  wave  for  M^-l,  and  everywhere  except  far  upstream 
for  M^>>1 . Its  actual  utility  can  only  be  judged  a posteriori. 
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Based  on  the  above  developments,  the  perturbation  expansion 
can  now  proceed  by  rewriting  (III-32)  and  (III-33)  for  i = 2 in 
terms  of  R , R^,  and  £,  defined  by  T,  - K2^l^*z  and  then  seek 
Rp  and  R^,  in  successive  approximation  by  writing 
fn\  m 9 to\ 

••  (III-34) 
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(III— 35) 


Since  k,,  is  brought  into  the  differential  equations  (III— 32 ) and 
(III— 33)  when  they  are  written  in  terms  of  £,  the  c expansion  for 
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In  zeroth-order  (II 1—32 ) and(III-33)  are  satisfied  identi- 
cally. In  first-order,  both  (III-32)  and  III-33)  yield  the  same 
relationship 

Rp0)(c)  - *40)(O. 

In  second-order  (IiI-32)  and  ( 1 11—33)  yield,  respectively, 
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and 
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For  clarity,  the  argument  £ will  be  dropped  throughout 
the  remainder  of  this  report . 


Using  (III-37)  to  express  R,^  in  terms  of  , the  last 
two  equations  may  be  looked  upon  as  involving  R,^  , R^  , and 


R 


(0) 


The  difference  between  (III-38)  and(III-39)  is  a first- 

,(0) 


order  differential  equation  for  R^  ' in  the  form 
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(111-40) 


Forming  the  sum  of  (111-38)  and  (III-39) , the  derivative  of 
Rp^  may  be  eliminated  using  (III-40)  to  give  an  algebraic  equation 
for  R^1^  as 
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(III-41) 


For  convenience  in  calculating  the  asymmetry  of  the  normal- 
ized density  profiles,  the  origin  of  the  coordinate  system  is  taken 
where  (C.)  = 0.5.  With  this  specification,  the  integration  of 


the  second-order  result  gives 

1 + Tanh 
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'P 
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and  from  (I 11-37) 
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(III-42) 
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Equation  (111-41)  will  yield  immediate  results  for  once 

has  been  determined  from  the  succeeding  perturbation  calculations. 

The  third-order  expressions  obtained  from  (III— 32)  and 
(111-33)  introduce  the  next  higher  order  of  the  perturbation 


quantities.  The  results  are 
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Eliminating  R^^  and  by  means  of  (III— 37 ) and  (III-41), 

and  taking  the  difference  between  (111-44)  and  (IJ.I-45) , one  finds 


+ [2R<°>-l]  R«>  Jr<0>|2- 

i P ) p bn)*  L p . 


_61  R(0) 
4 (7)2  L P 


( I II— 46) 


An  expression  for  R^  ' is  found  using  (III-46)  in  the  sum 


of  (III-44)  and  (III-45).  The  result  is 
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Proceeding  as  before,  the  fourth-order  perturbation 
calculations  give 
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In  fifth-order,  the  results  are 
, (3) 
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and 
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There  is  a pattern  in  the  preceding  calculations,  although  it 
will  probably  have  escaped  the  reader's  attention.  As  is  pointed 
out  in  Appendix  C,  not  all  of  the  numerical  coefficients  in  the 
differential  equations  for  R^  are  independent,  nor  are  all  the 
coefficients  in  the  algebraic  results  for  R^^  independent.  The 
pattern  involves  a limited  number  of  numerical  coefficients  denoted 


by  A . 
n 
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Moreover,  as  also  pointed  out  in  Appendix  C,  there  is  a 
conveniently  economical  way  to  determine  the  as  functionals 
of  (1  - Rp0)).  One  finds 


(111-54) 

In  ( I 11—53)  and  (III-54),  the  coefficients  B and  C are 

n n 


functions  of  the  generalized  coefficients  and  the  integration 
constants  Q..  The  expressions  for  the  B and  C in  terms  of  the 
A^  are  also  given  in  Appendix  C.  Once  R^  has  been  determined 
for  each  value  of  £ via  ( I XI— 4 2) , the  corresponding  perturbed 
quantities  and  R^^  for  j = 1,  2,  3 may  be  evaluated  to 

determine  the  individual  contributions  to  the  normalized  density 


and  temperature  for  each  of  the  above  orders  of  the  perturbation 
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development.  The  density  and  temperature  profiles  depend  on 
the  Mach  number  through  the  R^,  R,p,  and  K expansions  in  terms  of 
the  Mach  number  parameter  e. 

The  discussion  of  the  Navier-Stokes  perturbation  results 
and  comparison  with  Runge-Kutta  numerical  integration  results  is 


deferred  until  Chapter  V. 


CHAPTER  IV 


METHODS  TO  SOLVE  THE  HIGHER  ORDER 
GASDYNAMIC  EQUATIONS 

Numerical  Integration  of  the  Higher 
Order  Equations 


As  indicated  earlier  in  this  report  serious  ditnculties 

arise  when  attempts  are  made  numerically  to  integrate  the  Burnett, 

Super-Burnett,  or  Super-Super-Burnett  equations.  The  nature  of 

the  difficulties  may  be  understood  from  an  analysis  of  the  singular 

points  of  the  differential  equations  using  the  linearized  equations, 

just  as  in  the  preceding  chapter.  Following  the  Navier-Stokes 

development,  the  variations  of  the  normalized  flow  velocity  and 

normalized  temperature  from  their  Rankine-Hugoniot  values  are 

<iz 

assumed  to  be  proportional  to  e far  upstrea-n  and  far  downstream. 
The  coupled  pairs  of  linear  differential  equations  then  yield  char- 
acteristic polynomials  for  the  k^. 

The  (linearized)  Burnett  equations  yield  a fourth  degree 
characteristic  polynomial. 


4 3 2 2 
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= 0, 


(IV-1) 
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(IV— 2 ) 
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The  roots  of  (IV-1)  are  presented  graphically  as  a function  of 
in  Appendix  B. 

The  upstream  roots  indicate  that  numerical  integration  is 
feasible  starting  downstream,  as  in  the  Navier-Stokes  case,  provided 
M^<1.89.  For  between  1.0  and  approximately  1.4,  the  roots  of 
the  upstream  characteristic  polynomial  are  real  and  positive  or  are 
complex  conjugates  with  positive  real  parts.  The  solutions  of  the 
linearized  equations  are  thus  damped  or  oscillatory  and  damped,  and 
the  upstream  Rankine-Hugoniot  asymptotic  values  are  attained  for  all 
solutions  as  z ->  - <»,  For  values  of  between  approximately  1,4 
and  approximately  1.9,  all  roots  are  complex  with  positive  real 
parts.  All  solutions  thus  exhibit  a somewhat  stronger  oscillatory 
behavior  far  upstream  for  between  approximately  1.4  and  approxi- 
mately 1.9,  but  remain  damped  as  z -+  - “.  Above  = 1.9,  the  real 
part  of  one  pair  of  conjugate  roots  changes  sign,  and  there  is  a 
manifold  of  undamped  solutions  which  precludes  numerical  integration 
starting  downstream. 

The  roots  of  the  downstream  Burnett  characteristic  polynomial 
preclude  numerical  integration  starting  upstream  for  all  M^. 
Therefore,  numerical  integration  proceeding  in  the  upstream  direc- 
tion is  the  only  feasible  choice.  However,  the  integration  may  be 
expected  to  be  successful  for  only  the  range  of  Mach  numbers  indi- 
cated above,  and  experience  confirms  this  expectation. 

The  (linearized)  Super-Burnett  equations  yield  a sixth 
degree  characteristic  polynomial, 

;< 
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with  defined  as  before. 

The  results  for  the  upstream  and  downstream  roots  as  func- 
tions of  are  given  in  Appendix  B.  This  graphical  presentation 
shows  that  some  of  the  real  roots  and  some  of  the  real  parts  of 
the  complex  roots  have  opposite  sign,  both  upstream  and  downstream, 
for  all  Mach  numbers.  Under  these  circumstances  there  is  no  reason 
to  expect  numerical  integration  to  succeed  proceeding  in  either 
direction  for  any  Mach  number,  and  experience  confirms  this  expec- 
tation 

The  results  of  the  linearization  of  the  Super-Super- 
Burnett  equations  are  qualitatively  the  same  as  those  found  for 
the  Super-Burnett  equations.  The  are  determined  from  the 
following  characteristic  polynomial : 
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fe 
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The  graphical  representation  of  the  roots  of  this  characteristic 
polynomial,  given  in  Appendix  B,  reveals  that  some  of  the  real 
roots  and  some  of  the  real  parts  of  the  complex  roots  again  have 
opposite  sign,  both  upstream  and  downstream,  for  every  Mach  number. 

Since  the  conventional  numerical  methods  fail  when  applied 
to  the  higher  order  equations,  the  remainder  of  this  chapter  is 
devoted  to  the  perturbation  method  developed  in  the  previous  chapter. 


Perturbation  Expansion  Applied  to  the 
Burnett  Equations 


The  burnett  equations  containing  + P^  and  q^  + 

(2) 

q , may  be  analyzed  in  the  same  way  as  the  Navier-Stokes  equations 
to  determine  a suitab] e Mach  number  dependent  length  scale.  The 
results  of  such  an  analysis  are  qualitatively  the  same  as  for  the 
Navier-Stokes  case,  and  thus  suggest  a perturbation  expansion  which 
incorporates  the  Burnett  result  for  the  physically  significant 
(negative,  tending  to  zero  as  M^^00).  The  expansion  of 
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Since  the  first  term  in  (IV-5)  is  the  same  as  the  first 
term  in  the  Navier-Stokes  ^ expansion  (111-36),  the  Burnett  calcu- 
lations do  not  effect  an;  changes  to  the  Navier-Stokes  perturbation 
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results  through  second  order  in  £,  although  of  course  the  numerical 
values  of  k are  different. 

Following  the  same  procedures  used  for  the  Navier-Stokes 
equations,  the  Burnett  calculations  in  third  order  yield 
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Referring  to  equations  (111-46)  and  (III-47),  it  can  be  seen  that 
the  third-order  perturbation  results  for  the  Navier-Stokes  equations 
and  the  Burnett  equations  are  of  the  same  form;  only  the  values  of 
the  corresponding  coefficients  have  changed. 

The  fourth-order  Burnett  calculations  yield 
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Again,  the  form  of  the  Burnett  results  is  the  same  as  those  found  in 
fourth  order  for  Navier-Stokes,  except  that  the  numerical  coeffi- 
cients are  different. 

The  fact  that  the  perturbation  equations  have  the  same  form 
for  the  Burnett  equations  as  the  Navier-Stokes  equations  applies 
also  to  the  Super-Burnett  equations  and  the  Super-Super-Burnett 
equations.  This  "pattern"  is  what  gives  general  utility  to  the 
material  in  Appendix  C concerning  integration  of  th<  perturbation 
equations. 

The  fifth-order  Burnett  results  are 
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Discussion  of  the  solutions  to  the  perturbation  equations 
is  deferred  until  Chapter  V. 


Perturbation  Expansion  Applied  to  the 
Super-Burnett  Equations  and  the 
Super-Super-Burnett  Equations 


The  Super-Burnett  equations,  containing  + P^  + P^ 
and  q^  + q^  4-  q^,  constitute  two  coupled  third-order  non- 
linear differential  equations  for  the  velocity  (or  density)  and 
temperature.  The  equations  provide  a rationale  for  a perturbation 
expansion  which  incorporates  the  Super-Burnett  result  for  the 
physically  significant  (negative,  tending  to  zero  as  M^-*00), 

just  as  the  Navier-Stokes  and  Burnett  equations  provide  such  a 

^2^2  2 

rationale.  The  expansion  for  — — <2  in  powers  of  e = Mj  - 1 


follows  from  the  characteristic  polynomial  (IV-4)  , 
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Since  the  first  two  terms  of  (IV-12)  are  the  same  as  the  corres- 
ponding i wo  terms  in  the  Burnett  expansion,  the  Super-Burnett 
calculations  do  not  change  the  Burnett  perturbation  results  thro >gh 
third  order  in  e,  except  of  course  through  the  numerical  value  of 

K2  ’ 

In  fourth  and  fifth  orders  in  e,  the  calculations  yield 
equations  for  the  normalized  density  and  temperature  which  are  of 
the  same  form  as  the  perturbation  equations  derived  in  the  Navier- 
Stokes  and  Burnett  perturbation  work.  The  generalized  equations 
and  the  coefficients  appropriate  to  Super-Burnett  are  given  in 


Appendix  C. 
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The  length  scale  expansion  for  the  Super-Super-Burnett 


equations  is 


P2a2  f 3 

' ^2M2  L 7 


, 102  .2  14190  _3  , 1058952 

e + e e + ~ 

2 (7;  4 (7;  8(7) 


f s4  + 


(IV-13) 


It  is  clear  from  this  form  for  k ^ that  the  Super-Super-Burnett 
calculations  do  not  change  the  Super-Burnett  results  through  fourth 
order  in  e.  In  fifth  order,  the  results  are  of  the  same  form  as  the 
fifth-order  results  for  Navier-Stokes,  Burnett,  and  Super-Burnett. 
The  Super-Super-Burnett  coefficients  are  also  given  in  Appendix  C. 

The  solutions  of  the  perturbation  equations  derived  in 
this  section  for  me  normalized  density  and  temperature  profiles 
will  be  discussed  in  the  next  chapter. 

As  a concluding  remark  on  the  work  required  in  the  pertur- 
bation development,  the  determination  of  the  eighty-four  coeffi- 
cients needed  through  fifth  order  in  e was  a major  effort.  The 
pattern,  which  ultimately  evolved,  was  not  obvious  or  confirmed 
until  the  Super-Burnett  development  was  completed. 

The  intuitive  expectation  that  the  R^^  and  R,^ 

(j  = 1,2,3,»»*»)  should  vanish  when  = 1,  in  order  to  satisfy 

the  Rankine-Hugoniot  relations,  piovided  a convenient  mechanism 
for  checking  the  coefficients  of  the  powers  of  in  each  per- 

turbation equation.  That  is,  the  sum  of  the  coefficients  of  the 
powers  of  R^^  had  to  vanish  for  each  perturbation  equation  in 
order  to  satisfy  the  Rankine-Hugoniot  conditions  far  upstream. 
Although  this  was  a convenient  check  on  the  calculations,  the 
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mechanism  was  not  foolproof.  It  turns  out  that  within  the  arith- 
metic processes  of  substitution,  expansion  and  collecting  terms, 
it  is  possible  to  carry  through  a mathematical  error,  undetected; 
yet,  with  a result  that  satisfies  the  vanishing  of  this  coefficient 
sum.  The  calculations  are  lengthy  and  require  meticulous  attention 
to  preclude  introducing  such  casual  errors.  Errors  of  this  sort 
may  go  undetected  except  by  separate  and  independent  checks  of  the 
calculations. 

In  the  present  work,  alternate  approaches  were  employed  to 
determine  firm  values  of  the  coefficient  sets  through  A^,  A^ 
through  Ag,  A^  through  A12,  A^g  through  A^g,  and  A17  through  A2l 
as  these  are  the  coefficient  sets  of  the  powers  of  in  each 

of  the  gasdynamic  developments. 

The  final  values  for  the  coefficients  in  each  of  the  sets 


are  based  on  the  fact  that  the  same  results  were  obtained  from  at 
least  two  independent  calculations  for  each  order  of  e of  the 
perturbation  expansion  in  each  of  the  gasdynamic  developments. 


CHAPTER  V 


SHOCK  WAVE  PROFILES  BASED  ON  THE  METHODS  OF 
NUMERICAL  INTEGRATION  AND 
PERTURBATION  EXPANSION 


Results  for  the  Navier-Stokes  Equations  and 
the  Burnett  Equations  Based  on  the 
Method  of  Numerical  Integration 


An  increment  size  of  0.05  upstream  Maxwell  mean  free  paths 
was  chosen  for  the ’ fourth-order  Runge-Kutta  numerical  integration. 
Smaller  increment  'sizes  gave  the  same  results  for  normalized  density 
and  temperature  everywhere  within  the  shock  wave  (through  six 
decimal  places) . The  linearized  Navier-Stokes  equations  and  the 
linearized  Burnett  equations  were  used  to  determine  the  initial 
values  for  the  numerical  integration,  starting  from  the  assigned 
value  downstream 


= - 0.001  . 


Figures  3 through  6 exhibit  the  normalized  density  profiles 
and  normalized  temperature  profiles  for  M^  - 1.5  and  M^  = 2.0. 
Additional  results  for  the  Navier-Stokes  normalized  density  profile 
at  higher  Mach  numbers  are  given  in  the  next  section.  The  Navier- 
Stokes  equations  were  integrated  numerically  by  other  authors  many 


J 


Figure  3 - Gas  density  profiles  in  a normal  shock  wave  (M^  = 1.5) 
based  on  the  numerical  integration  of  the  Navier-Stokes 
and  Burnett  gasdynamic  differential  equations. 

Navier-Stokes  normalized  density 

Burnett  normalized  density 


Maxwell  Mean  Free  Paths 


Figure  4 - Gas  density  profiles  in  a normal  shock  wave  (M^  = 2.0) 
based  on  the  numerical  integration  of  the  Navier-Stokes 
and  Burnett  gasdynamic  differential  equations. 

Navier-Stokes  normalized  density 

Burnett  normalized  density 
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Figure  5 - Gas  temperature  profiles  in  a normal  shock  wave  (M  =1.5) 
based  on  the  numerical  integration  of  the  Navier-Stokes 
and  Burnett  gasdynamic  differential  equations. 

__  Navier-Stokes  normalized  temperature 


Burnett  normalized  temperature 


Z/A  - Maxwell  Mean  Free  Paths 
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Figure  6 - Gas  temperature  profiles  in  a normal  shock  wave  (M  =2.0) 
based  on  the  numerical  integration  of  the  Navier-Stokes 
and  Burnett  gasdynamic  differential  equations. 

Navier-Stokes  normalized  temperature 

''  Burnett  normalized  temperature 


Maxwell  Mean  Free  Paths 
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years  ago,  but  (as  mentioned  in  Chapter  I)  the  Burnett  equations 
have  not  been  integrated  numerically  before. 

Figure  7 shows  the  Mach  number  dependence  of  the  distance 
from  the  center  of  the  normalized  temperature  profile  to  the  center 
of  the  normalized  density  profile  (the  former  is  always  upstream 
relative  to  the  latter). 

As  Figures  3 through  7 indicate,  there  are  definite  differ- 
ences between  the  Navier-Stokes  shock  structure  and  the  Burnett 
shock  structure,  even  for  M^<2;  the  difference  in  Mach  number 
dependence  of  the  separation  of  the  profiles  is  particularly 
conspicuous.  The  differences  in  the  profiles  themselves  are  accen- 
tuated in  plots  of  the  asymmetry  quotients  (introduced  in  Chapter  I) 
as  functions  of  Mach  number,  as  shown  in  Figure  8 and  Figure  9. 

The  results  of  this  section  suggest  that  the  attention  of 
experimentalists  should  be  directed  toward  asymmetry  quotients  for 

/ 

,,  temperature  profiles  as  well  as  density  profiles,  and  toward  the  / 

separation  of  the  profiles. 

Results  for  the  Navier-Stokes  Density  Profile 
Based  on  the  Methods  of  Numerical 
Integration  and  Perturbation 
Expansion 

The  success  of  the  perturbation  method,  as  applied  to  the 
Navier-Stokes  density  profile,  may  be  judged  from  the  comparisons 
in  Figures  10  through  13.  It  is  clear  from  the  figures  that  the 


Maxwell  Mean  Free  Paths 
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Figure  7 - Mach  number  (M^)  dependence  of  the  temperature-density 

separation  in  a normal  shock  wave  based  on  the  numerical 
integration  of  the  Navier-Stokes  and  Burnett  gasdynamic 
differential  equations. 

Az  _ z z 

Profile  separation  — y = y ~ y 

1 *1  Rp  = 1/2  1 Rt  = 1/2 

Navier-Stokes 
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Burnett 


Density  Asymmetry  Quotient 
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Figure  1U  - Navier-Stokes  density  profiles  in  a normal  shock  wave 
(M^  = 1.5).  A comparison  of  density  profiles  based  on 
the  methods  of  numerical  integration  and  perturbation 
expansion. 

Numerical  integration 

Zeroth  order  perturbation 

First  order  perturbation 

Second  order  perturbation 
Third  order  perturbation 

The  third  order  perturbation  profile  is  indiscernible 
from  the  numerical  integration  profile  within  the  scale 
of  the  figure. 
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Figure  11  - Navier-Stokes  density  profiles  in  a normal  shock  wave 
(M^  =2.0).  A comparison  of  density  profiles  based  on 
the  methods  of  numerical  integration  and  perturbation 
expansion. 
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Zeroth  order  perturbation 

First  order  perturbation 
— Second  order  perturbation 


Third  order  perturbation 
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Figure  12  - Navier-Stokes  density  profiles  in  a normal  shock  wave 
(M^  = 3.0).  A comparison  of  density  profiles  based  on 
the  methods  of  numerical  integration  and  perturbation 
expansion. 


Numerical  integrati”" 

Zeroth  order  perturbation 

First  order  perturbation 
Second  order  perturbation 


Third  order  perturbation 
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Figure  13  - Navier-Stokes  density  profiles  in  a normal  shock  wave 
(M^  =4.0).  A comparison  of  density  profiles  based  on 
the  methods  of  numerical  integration  and  perturbation 
expansion. 


Numerical  integration 
Zeroth  order  perturbation 
First  order  perturbation 
Second  order  perturbation 


Third  order  perturbation 


Mean  Free  Pa 


successive  approximations  of  the  perturbation  method  are  converging 
toward  the  Runge-Kutta  result.  The  agreement  between  the  pertur- 
bation results  and  numerical  integration  is  satisfactory  for  the 
lower  Mach  numbers  throughout  the  profile,  but  it  would  be  desirable 
to  know  the  perturbation  results  to  another  order  or  two  for  the 
higher  Mach  numbers  to  do  justice  to  the  upstream  portion  of  the 
profile. 

The  convergence  of  the  perturbation  results  toward  the  Runge- 
Kutta  result  is  also  manifested  in  plots  of  the  asymmetry  quotient, 
as  shown  in  Figure  14. 

Results  for  the  Burnett  Density  Profile  Based 
on  the  Methods  of  Numerical  Integration 
and  Perturbation  Expansion 

The  results  of  the  perturbation  method,  as  applied  to  the 
Burnett  density  profile,  are  compared  with  the  Runge-Kutta  results 
in  Figures  15  and  16.  Figures  17  and  18  give  the  perturbation 
results  for  Mach  numbers  beyond  the  range  of  applicability  of 
numerical  integration.  Figure  19  contains  a comparison  between  the 
asymmetry  quotient  as  determined  by  numerical  integration  (up  to 
M^=  2),  and  the  assymmetry  quotient  as  determined  by  perturbation 

4. 

expansion. ' 

Figures  15  and  16  show,  as  in  the  case  of  Navier-Stokes 

equations,  that  the  perturbation  results  are  converging  toward  the 

Z /Ay 

'since  Rp  is  symmetric  yielding  Qp=  1 for  all  Mach 
numbers,  these  results  have  been  omitted  from  all  plots  of  the 
asymmetry  quotient. 


I 
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1 gure  14  - Density  asymmetry  quotient  (Q  ) versus  M from  the  Navier- 
Stokes  perturbation  expansion.  A comparison  of  the  results 
for  Q based  on  the  methods  of  perturbation  expansion  and 
numerical  integration. 

Numerical  integration 


First  order  perturbation  expansion 
Second  order  perturbation  expansion 
Third  order  perturbation  expansion 


Figure  15  - Burnett  density  profiles  in  a normal  shock  wave  (M^  = 1.5). 

A comparison  of  density  profiles  based  on  the  methods  of 
numerical  integration  and  perturbation  expansion. 

Numerical  integration 

First  order  perturbation 

Second  order  perturbation 

_ Third  order  perturbation 

The  third  order  perturbation  profile  is  indiscernible  from 
the  numerical  integration  profile  within  the  scale  of  the 
figure . 


Figure  16  - Burnett  density  profiles  in  a normal  shock  wave  (M^  = 2.0). 

A comparison  of  density  profiles  based  on  the  methods  of 
numerical  integration  and  perturbation  expansion. 

Numerical  integration 

First  order  perturbation 

Second  order  perturbation 

_ Third  order  perturbation 
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Figure  17  - Burnett  density  profiles  in  a normal  shock  wave  (M.  = 3.0). 
A comparison  of  density  profiles  determined  from  tne 
perturbation  expansion. 

__  First  order  perturbation 

Second  order  perturbation 

___  Third  order  perturbation 


Mean  Free  Paths 
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Figure  Id  - Burnett  density  profiles  in  a normal  shock  wave  (M  = 4.0). 
A comparison  of  density  profiles  determined  from  tne 
perturbation  expansion. 

First  order  perturbation 

. Second  order  perturbation 

Third  order  perturbation 
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1.0  2.0  3.0  4.0  M1 

Figure  19  - Density  asymmetry  quotient  (Q  ) versus  from  the 

Burnett  perturbation  expansioR.  A comparison  of  the 
results  for  Q based  on  the  methods  of  perturbation 
expansion  and^numerical  integration. 

Numerical  integration 

First  order  perturbation  expansion 

# Second  order  perturbation  expansion 

Third  order  perturbation  expansion 


Runge-Kutta  result.  This  conclusion  is  tempered,  but  not  vitiated, 
by  the  comparison  of  asymmetry  quotients  in  Figure  19,  which 
likewise  indicate  a tendency  toward  convergence. 

Intuitively  Trustworthy  Parts  of  the  Results 
for  the  Density  Profile  through  Super- 
Super-Burnett  Order  Based  on  the 
Method  of  Perturbation  Expansion 

The  Chapman-Enskog  development  provides  a criterion  for 
identifying  the  trustworthy  parts  of  predictions  of  the  Navier- 
Stokes  or  higher  order  gasdynamic  equations.  Predictions  of  the 
Navier-Stokes  equations  may  be  trusted  only  through  first  order  in 
the  mean  free  path,  because  the  Burnett  contributions  to  the  stress 
and  heat  flux  are  of  second  order  in  the  mean  free  path.  Similarly, 
predictions  of  the  Burnett  equations  may  be  trusted  only  through 
second  order  in  the  mean  free  path,  because  the  Super-Burnett  con- 
tributions to  the  stress  and  heat  flux  are  of  third  order  in  the 
mean  free  path,  etc. 

In  the  previous  chapter  it  was  pointed  out  that  the  pertur- 
bation expansion  may  be  applied  with  equal  justification  to  all  the 
higher  order  gasdynamic  equations  through  Super-Super-Burnett  order, 
and  presumably  this  feature  persists  to  all  orders  of  the  Chapman- 
Enskog  development.  Accordingly,  the  results  of  a particular  order 
of  the  perturbation  expansion  may  be  regarded  as  trustworthy  only 
if  they  are  unaffected  by  additional  Chapman-Enskog  contributions 
to  the  stress  and  heat  flux.  Thus,  the  Navier-Stokes  result  for 

R^  is  trustworthy,  but  R^\  R^  and  R^  are  not  trustworthy 

P P P p 
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because  they  are  modified  by  the  inclusion  of  Burnett  (and  higher 

Chapman-Enskog  order)  contributions  to  the  stress  and  heat  flux. 

Likewise,  the  Burnett  results  for  R^  and  R^  are  trustworthy, 

(2)  (3) 

but  Rp  and  Rp  are  not  trustworthy,  etc. 

According  to  the  viewpoint  explained  in  the  preceding  two 
paragraphs,  the  trustworthy  part  of  the  Super-Burnett  result  for 
the  density  profile  is 

R<°>  f E R<»  + E2  R«) 

P P P 

and  the  trustworthy  part  of  the  Super-Super-Burnett  result  is 

R<°>  + c R<3>  + e2R<2>  + e3R<3>  . 

P P P P 

These  trustworthy  parts  of  the  predictions  of  the  density  profile 
are  exhibited  in  Figures  20  through  23,  and  the  corresponding  results 
for  the  asymmetry  quotient  are  exhibited  in  Figure  24.  The  trust- 
worthy parts  of  the  Navier-Stokes  and  Burnett  calculations  have  been 
omitted  from  Figures  20  through  24  to  avoid  undue  cluttering  of  the 
figures,  but  the  Navier-Stokes  numerical  integration  results  have 
been  included  to  provide  an  indication  of  the  magnitude  of  higher 
order  Chapman-Enskog  effects.  These  figures  may  be  said  to  provide 
a hint  that  extremely  high  orders  of  uhe  Chapman-Enskog  development 
will  not  be  required  for  a satisfactory  theory  of  shock  structure. 


Figure  20  - Super-Burnett  and  Super-Super-Burnett  density  profiles 
in  a normal  shock  wave  (M^  = 1.5).  A comparison  of  the 
trustworthy  parts  of  the  perturbation  results  for  Super- 
Burnett  and  Super-Super-Burnett  with  the  Navier-Stokes 
numerical  integration  results. 

__  Navier-Stokes  numerical  integration 

, , . Second  order  Super-Burnett  perturbation 

Third  order  Super-Super-Burnett  perturbation 


Z/A  - Maxwell  Mean  Free  Paths 


Figure  21  - Super-Burnett  and  Super-Super-Burnett  density  profiles 
in  a normal  shock  wave  (M^  = 2.0).  A comparison  of  the 
trustworthy  parts  of  the  perturbation  results  for  Super- 
Burnett  and  Super-Super-Burnett  with  the  Navier-Stokes 
numerical  integration  results. 

' Navier-Stokes  numerical  integration 

# # Second  order  Super-Burnett  perturbation 
Third  order  Super-Super-Burnett  perturbation 


Z/A,  - Maxwell  Mean  Free  Paths 
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Figure  22  - Super-Burnett  and  Super-Super -Burnett  density  profiles 
in  a normal  shock  wave  (M^  =3.0).  A comparison  of  the 
trustworthy  parts  of  the  perturbation  results  for  Super- 
Burnett  and  Super-Super-Burnett  with  the  Navier-Stokes 
numerical  integration  results. 

Navier-Stokes  numerical  integration 

, . Second  order  Super-Burnett 

Third  order  Super-Super-Burnett 


Figure  23  - Super-Burnett  and  Super-Super-Burnett  density  profiles 
in  a normal  shock  wave  (M^  =4.0).  A comparison  of  the 
trustworthy  parts  of  the  perturbation  results  for  Super- 
Burnett  and  Super-Super-Burnett  with  the  Navier-Stokes 
numerical  integration  results. 

Navier-Stokes  numerical  integration 

. .Second  order  Super-Burnett 

Third  order  Super-Super-Burnett 


Figure  24  - Density  asymmetry  quotient  (Q  ) versus  M1  from  the  Super- 
Burnett  and  Super-Super-Bu^ne^t  perturbation  expansions. 

A comparison  of  the  trustworthy  parts  of  the  perturbation 
results  for  Super-Burnett  and  Super-Super-Burnett  expan- 
sions with  Navier-Stokes  numerical  integration  results. 

Navier-Stokes  numerical  integration 

, Second  order  Super-Burnett  perturbation 

_ Third  order  Super-Super-Burnett  perturbation 


CHAPTER  VI 


COMMENTS  AND  OUTLOOK 


The  general  pattern  among  the  numerical  coefficients  in  the 
perturbation  results  derived  in  the  preceding  chapters  may  be 
obtained  from  the  differential  equation 
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and  the  algebraic  relation 
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by  expanding  R„  and  R„  in  powers  of  c.  This  observation  constitutes 
P l 

an  explanation  or  interpretation  of  the  general  pattern.  The 
first-order  differential  equation  (VI-1)  is  in  terms  of  the  unex- 
panded normalized  density  R^  and  may  be  integrated  using  standard 
Runge-Kutta  numerical  integration  methods.  The  normalized 
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temperature  profile  may  then  be  determined  from  the  algebraic 
equation  (VI-2)  relating  and  R^. 

In  carrying  out  the  algebraic  determination  of  the  R^' 
for  an  expansion  of  k^.  in  which  the  numerical  values  of  the  coeffi- 
cients were  unspecified 


(j) 


_ P2a2 
<2  = *2M2 


i 2 3 4 

(o^e  + a2e  + a3e  + a^e  + 


(VI-3) 


it  was  found  that  all  of  the  terms  containing  the  expansion  coef- 
ficients (a^)  cancel  during  the  calculation.  Thus,  the  final  rela- 
tionship (VI-2)  between  R^  and  R^is  independent  of  the  choice  for 
This  is  a satisfying  result  because  R^  should  be  uniquely 
determined  by  R^,  assuming  each  is  a monotonic  function  of  z, 
quite  independent  of  any  reciprocal  length  scale  k,,.  Also,  in  the 

analytic  expression  for  the  R^\  it  turns  out  that  the  a.  are 

p i 

uniquely  determined  for  each  set  of  the  gasdynamic  equations  from 

the  requirement  that  the  ratio  R^/R^  be  bounded  throughout  the 

P P 

shock  wave  for  each  (j).  In  order  for  this  ratio  to  be  finite,  the 
must  have  the  values  as  specified  in  Table  I for  each  of  the 
gasdynamic  developments. 

The  upstream-downstream  symmetry  observed  in  the  basic 

gasdynamic  differential  equations  leads  to  interrelationships 

among  the  perturbation  coefficients  A in  (VI-1).  Since  formally 

n 

2 

equivalent  expansions  for  dR  /dz  may  be  obtained  based  on  (M  -1) 

P 2 

2 

and  <2  or  (M“-l)  and  k^,  and  since  the  expansions  must  represent 
one  and  the  same  differential  equation  for  R , roughly  half  the 
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TABLE  I 

COEFFICIENTS  « FOR  <2 
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coefficients  in  each  order  of  £ are  determined  by  the  other  half 
(and  coefficients  from  lower  orders  of  e). 

The  consequences  of  upstream-downstream  symmetry  indicate 
that  it  may  be  possible  to  extend  the  perturbation  results  to  higher 
order  in  e with  somewhat  less  labor  than  was  expended  in  this 
report . 

All  profile  and  asymmetry  data  calculated  in  this  report. 

have  been  for  Maxwell  molecules,  for  which  the  viscosity-temperature 

relation  is  y^T.  Numerical  integration  results  for  the  Navier- 

Stokes  gasdynamic  equations  may  also  be  obtained  for  more  realistic 

intermolecular  potentials  by  using  the  viscosity-temperature 

relation  y«T  , where  w is  fractional.  For  purposes  of  these  final 

remarks,  the  Navier-Stokes  results  for  a realistic  potential 
0 68 

(y«T  ’ ) have  been  calculated  and  are  compared  with  recent  experi- 

mental data  and  perturbation  results  for  Maxwell  molecules  in 
Figure  25. 

The  perturbation  results  for  the  asymmetry  in  the  density 
profile  in  c Maxwell  gas,  although  based  on  an  idealized  molecular 
model, resemble  the  trend  of  discrepancy  between  current  experimental 
results  and  the  quantitatively  incorrect  Navier-Stokes  results. 

It  seems  reasonable  to  hope  that  use  of  a more  realistic 
potential  model,  together  with  an  extension  of  the  perturbation 
method  to  one  or  perhaps  two  additional  orders,  might  lead  to 
definitive  checks  on  the  nonlinear  features  of  the  Boltzmann 
c<  llision  operator,  especially  in  view  of  the  inadequacy  of  the 
Navier-Stokes  equations  even  for  low  Mach  numbers. 


Ill 


Figure  25  - Density  asymmetry  quotient  (Q  ) in  a normal  shock  wave 

as  a function  of  upstream  Macn  number  A comparison 

of  Navier-Stokes  numerical  integration  (realistic  poten- 
tial) and  third  order  perturbation  (Maxwell  molecules) 
results  with  experiment. 

. — . Navier-Stokes  numerical  integration 
_ Navier-Stokes  third  order  perturbation 

Super-Super-Burnett  third  order  perturbation 

O Experiment  - Alsmeyer  (1975) 

A Experiment  - Schmidt  (1969) 
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APPENDIX  A 


THE  MAXWELL  EIGENFUNCTIONS 


The  eigenfunctions  ^(c)  used  in  the  moment  calculations 

of  Chapter  II  are  formed  from  associated  Laguerre  polynomials 

.£+1/2,  2.  , , , £ „ . + 

L (c  ) and  zonal  harmonics  c P.(c  ), 
r 36  Z 


0ft\  _ TX+1/2  t 2\  1 T>  X 

r» £(c)  — (c  ) c P£(cz)  • 


(A-l) 

The  associated  Laguerre  polynomials  may  be  determined  from  the 
generating  function 
2 t_ 

- 00 

n 0J-1  /O  0 ^ 

(A- 2) 


-c 


1-t 


l L*+1/2  (c2) 


(l-t)£+3,/2  r=0  r 

and  the  Legendre  polynomials  P (x)  may  be  determined  from 
Rodrigues'  formula 


P (x)  = 


2nn ! dxn 


x2-l 


(A-3) 


The  Maxwell  eigenfunctions  needed  in  this  effort  together 

with  their  associated  eigenvalues,  c .,  are  given  in  Table  II. 

r , a 

In  Table  II,  A?  is  a pure  number  which  has  been  defined 
in  Chapter  II. 

-f 

Useful  relations  among  the  ¥ (c)  are  listed  below: 

r , 36 


cz  ’i'o,0(c)  = 


C**r,0<c>  -*rfl(c)  “ Vl,l<c>’  ^ 


(A— 4) 
(A-5) 


■fWe  adopt  the  convention  that  A denotes  a unit  vector 
parallel  to  A for  any  A. 


117 


l 

< 

15 

\ 


V 


i 

4 


b 


| 

si 


TABLE  II 

MAXWELL  EIGENFUNCTIONS  AND  ASSOCIATED  EIGENVALUES 
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'r  ,1 


I-'2 


15  5 2 1 4 

8 - 2 c + 2 c 
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\ APPENDIX  B 

J ROOTS  OF  THE  CHARACTERISTIC  POLYNOMIALS 
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The  Mach  number  dependent  length  scales  alluded  to  in 
Chapters  III  end  IV  (in  conjunction  with  the  linearization  of  the 
gasdynamic  differential  equations)  may  be  determined  by  solving 
for  the  roots  of  the  characteristic  polynomials  which  are  obtained 
as  a result  of  the  linearization  process.  The  characteristic 

^iMi 

polynomials  in  L.  = ■-  — «.  given  below,  containing  the  Mach 

Piai  1 

number  M^  as  a parameter,  were  so  obtained. 

NAVIER-STOKES: 


20  + (9-23M^)  Li  + 6 mJ  [mJ  - l)  - 0 


(B— 1) 


BURNETT: 


70  L 4 - 180  + 291  mJ 


81  - 207  M7 


m; 


Li  + 54  mJ  (mJ-1)  = 0 


(B-2) 


SUPER- BURNETT: 


7850  Lj  + 29160  - 28980  M,  L4 

i i i l 


20475  M4  - 5427  M^ 


- 36(291)  mJ  I? 


36(207)  - 36(81)  mJ 


- 54(36)  M, 


mJ-i 


(B—  3) 


SUPER-SUPER-BURNETT  : 


328399(625)L®  - 200(63) (521)lJ 


+ 6300(3791)M^  L*  - 40 (81) (49) (59)mJ  Ljj 


+350(27) (2287)M*  - 54(49)mJ 


.6  ,2 


2275Mi  - 603| 


+ 8(49)  (81)  (97)M^  - 8(49)  (243)M][  23M^  - 9 


+ 16(49) (729)M 


8 


mJ  - 11 


0 . 


(B-4) 


The  roots  of  each  characteristic  polynomial  far  upstream 
(L^)  or  far  downstream  (L«)  may  be  determined  for  any  selected 
upstream  Mach  number  (M^)  or  its  corresponding  Rankine-Hugoniot 
downstream  Mach  number  (M^)  using  standard  computer  library  routines. 
The  results  given  in  the  tables  and  figures  which  follow  were 
computed  to  twelve  decimal  accuracy. 

Figures  26,  27,  28,  and  29  present  plots  of  the  results  for 
versus  M^  and  versus  M^  as  determined  from  computer  solutions 
of  equations  (B-l) , (B— 2 ) , (B-3),  and  (B-4),  respectively.  The  real 
L_^  and  the  real  parts  of  the  complex  are  plotted  as  solid  lines 
and  the  imaginary  parts  of  the  complex  are  plotted  as  broken  lines 
in  the  figures.  The  corresponding  real  and  imaginary  parts  of  the 
conjugate  pairs  are  appropriately  identified  in  each  figure. 

The  physically  significant  roots  <2»  which  provide  appro- 
priate length  scales  for  each  of  the  gasdynamic  developments,  are 
presented  for  a range  of  upstream  Mach  numbers  in  Table  III 


Figure  26  - Roots  of  the  Navier-Stokes  characteristic  polynomial. 

Upstream  roots  (L,)  and  downstream  roots  (L„)  versus  M, 


Figure  27  - Roots  of  the  Burnett  Characteristic  polynomial. 

Upstream  roots  (L  ) and  downstream  roots  (L2  ) 


TABLE  III 


PHYSICALLY  SIGNIFICANT  <2 


Super-Super' 


M1 

Naiver-Stokes 

Burnett 

Super-Burnett 

Burnett 

1.1 

-0.07167703 

-0.06685100 

-0.06633523 

-0.06612403 

1.2 

-0.12285635 

-0.10793995 

-0.10594205 

-0.10447577 

1.3 

-0.16089605 

-0.13425755 

-0.13066329 

-0.12708902 

1.4 

-0.19005288 

-0.15164829 

-0.14673941 

-0.14077175 

1.5 

-0.21295431 

-0.16341863 

-0.15753600 

-0.14923111 

1.6 

-0.23130542 

-0.17153184 

-0.16496566 

-0.15452634 

1.  7 

-0.24625774 

-0.17720053 

-0.17017150 

-0.15784456 

1.8 

-0.25861475 

-0.18119856 

-0.17386660 

-0.15989620 

1.9 

-0.26895253 

-0.18403361 

-0.17651194 

-0.16112073 

2.0 

-0.27769376 

-0.18604662 

-0.17841433 

-0.16179769 

2.1 

-0.28515479 

-0.18747177 

-0.17978304 

-0.16210907 

2.2 

-0.29157651 

-0.18847010 

-0.18076370 

-0.16217567 

2.  3 

-0.29714518 

-0.18915860 

-0.18145944 

-0.16207889 

2.4 

-0.  30200677 

-0.18961924 

-0.18194438 

-0.16187431 

2.5 

-0.30627706 

-0.18991189 

-0.18227243 

-0.16160018 

2.6 

-0.  31004889 

-0.19008036 

-0. 13248322 

-0.16128305 

2.  7 

-0.  31339746 

-0.19015700 

-0.18260620 

-0.16094136 

2.8 

-0.31638419 

-0.19016584 

-0.18266339 

-0.16058798 

2.9 

-0.31905971 

-0. 19012490 

-0.18267137 

-0.16023179 

3.0 

-0. 32146603 

-0.19004769 

-0. 18264270 

-0.15987893 
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TABLE  III 

(continued; 

« 

M1 

Naiver-Stokes 

Burnett 

Super-Burnett 

Super-Super- 

Burnett 

3.1 

-0.32363828 

-0.18994443 

-0.18258696 

-0.15953354 

3.2 

-0. 32560601 

-0.18982285 

-0.18251144 

-0.15919835 

3.3 

-0.32739421 

-0.18968880 

-0.18242175 

-0.15887508 

3.4 

-0.32902418 

-0.18954675 

-0.18232219 

-0.15856475 

3.5 

-0.33051413 

-0.18940005 

-0.18221604 

-0.15826785 

3.6 

-0.33187970 

-0.18925125 

-0.18210585 

-0.15798450 

3.7 

-0.33313441 

-0.18910227 

-0.18199355 

-0.15771459 

3.8 

-0.33428998 

-0.18895455 

-0.18188062 

-0.75745783 

3.9 

-0.33535662 

-0.18880913 

-0.18176819 

-0.15721379 

4.0 

-0.33634327 

-0.18866679 

-0.18165710 

-0. 15698200 

APPENDIX  C 


GENERAL  FORM  OF  THE  EQUATIONS  FROM  THE  PERTURBATION  DEVELOPMENT 
AND  GENERAL  FORM  OF  THE  SOLUTIONS  TO  THE  DIFFERENTIAL  EQUATIONS 


FOR  THE  NORMALIZED  DENSITY 
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The  perturbation  methods  applied  in  each  of  the  gasdynamic 
developments  of  Chapter  II  produce  a set  of  first  order  linear 
differential  equations  for  the  perturbed  density  constributions 

R(0)  r(D  r(2)}  r(3),  etc.,  and  a set  of  algebraic  relationships 

P P P P 

for  the  perturbed  temperature  contributions 

etc.,  in  terms  of  the  perturbed  densities.  In  the  course  of  the 
perturbation  calculations,  a pattern  in  the  analytical  results 
recurred  in  each  of  the  gasdynamic  developments. 

Thus,  the  perturbation  equations  for  the  normalized  density 
and  temperature,  applicable  to  any  of  the  gasdynamic  developments, 
are  of  the  following  general  form: 


(C-A) 
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n(°)  = R(0)  . 

®T  0 ’ 


(C— 5) 


4“  ■ Rp1>+A9  (l-*p0>  8p0)  ; 


(C-6) 


42)  - »fsK’Kvs(r  Mr  r ; 

L J L J (C-7) 

43)  = Rp}+A  |i-2R(0)]r(2)+A  -R(1)+2A11R(0)R(1)+3A  o[r(0)]  R1 
^ P 9 P p 10  p 11  p p 12[  p J ( 

- 4r  ] Wr-uko)  ] Ml  Mi  4(c=s 

»i4>  ■ Rp4:i+A3  [1-2Ro°>)K03>  +A10Rp2>+2AllRp0)Rr+3A12  Rp0>  \ 
+ Ai3RD1)+2Ai4Ri0)Ri1)+3Ais  RJ0)  r(1)+4A1  1r(0)1  r(1) 

13  p 14  p p 15>-  p J p 16  p p 

» . 


(C-8) 

psl2 


p 3 ! o rp  r“l0“p  ,4nll*p  *p  'Jn12\ 


, A P(1)+9A  .(OMDiW 

+ A.,  ~R_  +2A1/R^  +3At  c 

13  p 14  p p 15L 


R<°>  R(1) 

P P 


2AR^V24a  +3A  R(0)[r(1)]  +A  R^+A  - R(0) 

9 p p 11  p 12  p p 17  p 18  p 

J J L J 


+ A [r<0>13+A  JrHVa  »<°>1  5 

A19  Kp  A2CRp  A2IRp 


(C— 9 ) 


The  coefficients  An  determined  in  each  of  the  gasdynamic 

developments  are  listed  in  Table  IV. 

The  solution  of  the  differential  equations  for  R^ , Rp^ » 

(3) 

Rp  , etc. , may  be  expressed  as  functionals  of  the  quantity 

(1-Rp  ^ ) • In  each  instance,  the  differential  equations  can  be  arranged 

in  the  form 


-r--6 — + 2R^0)-1  R(n)  = I(n). 
dC  P P 


(C-10) 


131 


I 


1 

r^ 

m 

00 

co 

o 

lO 

VO 

o 

»H 

U 

"O 

r^ 

o 

CO 

CO 

00 

vO 

rH 

vO 

rH 

QJ 

i 

* — 1 

CO 

<r 

OV 

rH 

vO 

00 

CX  H 

co 

CO 

av 

<3- 

Ov 

Ov 

CO 

MT 

0 4J 

CM 

CO 

rH 

MT 

vO 

CO 

1 

co  d) 

1 

•<r 

Ov 

CM 

i F 

CM 

i — i 

O 

CN 

M W 

VO 

CM 

vO 

O 

o d 

1 

CM 

CM 

rH 

a pq 

| 

d 

CO 

4-» 

m 

CO 

co 

O 

-vT 

<r 

o 

rH 

H 

00 

r^> 

o 

co 

rH 

rH 

co 

CO 

vO 

H 

0) 

1 

rH 

co 

vO 

O 

<r 

O 

00 

c 

co 

CO 

o\ 

00 

00 

CM 

co 

<r 

Pi 

CM 

CO 

o 

rH 

vO 

l 

d 

1 

CM 

*sT 

CM 

O 

CQ 

o 

O 

r"» 

1 

o 

vO 

co 

r- 

m 

H 

rH 

CO 

*<f 

rH 

25 

d) 

1 

| 

o 

cx 

M 

d 

H 

CO 

< 

D 

O' 

W 

25 

O 

hH 

H 

r^ 

CO 

CO 

in 

oo 

<r 

o 

<r 

rH 

< 

co 

a> 

CM 

CO 

CO 

GV 

ov 

co 

vD 

rH 

PQ 

u 

1 

rH 

oo 

vO 

MJ- 

CO 

rH 

CM 

0Q 

AJ 

vO 

vO 

o 

CN 

CO 

Ov 

00 

NT 

15 

d) 

m 

CM 

c*. 

vO 

CO 

CO 

rH 

1 

H 

c 

i 

rH 

1 

0- 

rH 

m 

rH 

> 

P 

Pi 

CTN 

00 

<1- 

VO 

M 

W 

d 

CO 

rH 

CN 

<r 

PP 

CO 

1 

CN 

CO 

rH 

W 

1 

►J 

W 

« 

< 

H 

H 

P 

O 

P 

c 

rH 

CO 

<r 

r— 1 

o 

O 

VO 

vO 

rH 

< 

<D 

vO 

o 

av 

ov 

CO 

CN 

CM 

VO 

CO 

CN 

in 

cn 

n 

o 

CN 

OV 

CO 

0 

O 

<r 

<r 

tH 

vO 

vO 

CO 

rH 

H 

4J 

CM 

CO 

rH 

CM 

CO 

CO 

1 

2: 

CO 

1 

1 

ov 

m 

CO 

00 

w 

1 

rH 

co 

CO 

r*. 

M 

H 

rH 

o 

0) 

1 

M 

•H 

P 

> 

P 

03 

w 

25 

o 

o 

vO 

vO 

vO 

vO 

co 

<r 

MT 

<r 

/•v 

/-v 

/■"N 

CN 

/*N 

r^s 

r^. 

v_/ 

V/ 

r>. 

v-^ 

CO 

00 

00 

CO 

<r 

V-/ 

VO 

vO 

vO 

CN 

CM 

CN 

CN 

vO 

‘<3- 

*H 

rH 

rH 

rH 

rH 

tH 

rH 

00 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

o 

rH 

CN 

CO 

<r 

in 

vO 

r- 

00 

av 

rH 

< 

c 

< 

< 

< 

<2 

<2 

<2 

<2 

<2 

133 


An  equivalent  expression  for  (C-10)  is 


] R(0)fl-R(0)' 

•PL  P . 


H0,]l  *rK0>]  ' 


(C-ll) 


a?K0) 


Rv  y r 

_e L jL  .(■»  L _(o) 

i-r<0>](  « r o J • 


(C-12) 


Using  the  chain  rule  to  obtain  the  functional  F^[l-R^] 


rom  I | 

dJF<n>[l-R<°4  I 

1-RP  j 1 0 'll  fl-R«»] 


(C-13) 


yields 


F(,,)  1 


-R(o)].-rdfi-R;o>i 

p J p 

' w 


K0)]>r] 


(C-14) 


With  F'n^(l-R^]  known,  both  sides  of  equation  (C-12)  are  exact 
differentials  and  may  be  integrated  to  yield  R^  with  the  inte- 
gration constant  determined  from  the  conditions  that  at  £ = 0, 

r(h)  = A <5no  where  6n0  is  the  Kronecker  delta. 

P 2 

For  exam}.  ,,  the  differential  equation  for  R^^  is 


+ 2 R<°>-1  RrL)  = A 1-R(0)  I 
P P IP 

L J _ W J w 


The  results  for  1^  and  1-R^  are 

P 


- A, 


{ -[  >-*r]2  hr' 


F(1)  1-R^°^  = - A 
P 


dfl-R(0)l 


] K0)] 


(C— 15) 


(C— 16) 


(C— 17 ) 
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Thus,  using  (C— 17)  in  (012)  and  integrating,  the  result  is 


r(D  = R(0) 

p p 


1_Ep0>ll  - Alta 

fl-R(0)' 

P 

+ Ql 

1 J 1 

1 

(C-18) 


where  is  an  integration  constant.  Evaluating  from  the 

and  R(1 
P P 


boundary  conditions  on  R^  and  R^  at  £ = 0 yields 


Qx  = j A1  In  2. 


(C-19) 


(2)  (3) 

The  differential  equations  for  Rp  and  R^  may  be  solved 
using  the  same  procedures  described  above.  The  solutions  are  given 
in  equations  ( III—  5 3)  and  (III— 54) . The  differential  equation  for 
Rp^  can  be  integrated  without  using  this  method.  The  result  for 


.(0) 


is  given  in  equation  (II I— 42 ) . 


,0) 


Since  the  differential  equations  for  the  R^  have  the  same 


general  form  for  each  gasdynamic  development,  the  solutions  for  the 

RpH^  may  also  be  expressed  in  a general  form.  The  coefficients 

of  the  integrated  results  are  functions  of  the  original  The 

results  for  and  R^  are  relatively  simple.  However,  the 

(2)  (3) 

coefficients  in  R^  and  R^  are  much  more  involved,  namely 


r(2)  _ R(0) 

P P 


1--R 


(0) 


1-R 


(0) 


+ B2  &n 


1-R 


(0* 


+B„ 


+BC 


l-R 


1-R 


ri 


i n 


Zn 


1-R 


(01 


+ B4£n 


1-R 


(0) 


1-R 


(0) 


Q2 


(C-29) 
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V and  VI  respectively. 
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TABLE  V 

COEFFICIENTS  B FOR 

n p 


n 

B 

n 

i 

_A4 

4 - Mi + 

2 

+ 

2A4  + A^ 

3 

A2  - 
A1 

2AlQl 

4 

-1/2AJ 

2 

5 

A1 

/ 
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TABLE  VI 

COEFFICIENTS  C FOR 

n P 

n Cn 

1 5A1A4  + A^  + 2A1A3  - A?-  3Ag  + (2A^  + Ag  + 4A4)Q1 
+ AjQ*  - A1Q2  + 2Q1Q2 

2 A6  + 2A?  + 3Ag  + (-A  - 2A4)01  + AjQ, 

3 -2A3  - 2AxA3  - 6A1A4  + (-3A?  f 2k^  + 4A4)QX  + 2AXQ3  - 2A]LQ2 

4 AxA3  + 2AxA4  + 1/2A3Q1 

5 3/2A3  - 2AlA3-  4A3A4-  3aJq3 

6 1/2A,3  + 1/2 A + {-2k2.  - 3A  )Q  - 5/2A  qJ  + Q3 

x.  o 1 41  111 

7 2A3  + 3AxA4  + 5aJq  - 3A 

8 -5/2 A3  + 3A*Q 

9 -1/6A3 

10  A3 

11  -A3 


